Convection_Diffusion_Temperature.cpp

/****************************************************************************
* Copyright (c) 2026, CEA
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
*****************************************************************************/
 
#include <Convection_Diffusion_Temperature.h>
#include <Navier_Stokes_std.h>
#include <Probleme_base.h>
#include <Fluide_reel_base.h>
#include <Discret_Thyd.h>
#include <Frontiere_dis_base.h>
#include <Param.h>
#include <Transport_Interfaces_base.h>
#include <TRUSTTrav.h>
#include <SFichier.h>
#include <Domaine_VF.h>
#include <Matrice_Morse.h>
#include <Champ_Uniforme.h>
#include <Perf_counters.h>
 
Implemente_instanciable(Convection_Diffusion_Temperature,"Convection_Diffusion_Temperature",Convection_Diffusion_Temperature_base);
// XD convection_diffusion_temperature eqn_base convection_diffusion_temperature INHERITS_BRACE Energy equation
// XD_CONT (temperature diffusion convection).
// XD attr penalisation_l2_ftd bloc_lecture penalisation_l2_ftd OPT to activate or not (the default is Direct Forcing
// XD_CONT method) the Penalized Direct Forcing method to impose the specified temperature on the solid-fluid interface.
 
Sortie& Convection_Diffusion_Temperature::printOn(Sortie& is) const { return Convection_Diffusion_Temperature_base::printOn(is); }
 
Entree& Convection_Diffusion_Temperature::readOn(Entree& is)
{
  assert(la_temperature);
  assert(le_fluide);
  Convection_Diffusion_Temperature_base::readOn(is);
  //Nom unite;
  //if (dimension+bidim_axi==2) unite="[W/m]";
  //else unite="[W]";
  Nom num=inconnue().le_nom(); // Plan for passive scalar equations
  num.suffix("temperature");
  Nom nom="Convection_chaleur";
  nom+=num;
  terme_convectif.set_fichier(nom);
  //terme_convectif.set_description((Nom)"Convective heat transfer rate=Integral(-rho*cp*T*u*ndS) "+unite);
  terme_convectif.set_description((Nom)"Convective heat transfer rate=Integral(-rho*cp*T*u*ndS) [W] if SI units used");
  nom="Diffusion_chaleur";
  nom+=num;
  terme_diffusif.set_fichier(nom);
  //terme_diffusif.set_description((Nom)"Conduction heat transfer rate=Integral(lambda*grad(T)*ndS) "+unite);
  terme_diffusif.set_description((Nom)"Conduction heat transfer rate=Integral(lambda*grad(T)*ndS) [W] if SI units used");
  solveur_masse->set_name_of_coefficient_temporel("rho_cp_comme_T");
  return is;
}
 
void Convection_Diffusion_Temperature::set_param(Param& param) const
{
  Convection_Diffusion_Temperature_base::set_param(param);
  param.ajouter_non_std("penalisation_L2_FTD",(this));
}
 
int Convection_Diffusion_Temperature::lire_motcle_non_standard(const Motcle& un_mot, Entree& is)
{
  if (un_mot=="penalisation_L2_FTD")
    {
      is_penalized = 1;
      eta = 1.e-12;
      Cerr << " Equation_base :: penalisation_L2_FTD" << finl;
      Cerr << "Equation_base: modele penalisation ibc = "<< choix_pena << finl;
      Nom mot, nom_equation;
      int dim;
      is >> mot;
      Motcle accolade_fermee="}", accolade_ouverte="{";
      if (mot == accolade_ouverte)
        {
          is >> mot;
          while (mot != accolade_fermee)
            {
              nom_equation = mot;
 
              Cerr << "penalisation_L2_FTD equation: "<< nom_equation << finl;
              const Transport_Interfaces_base& eq = ref_cast(Transport_Interfaces_base,probleme().get_equation_by_name(nom_equation));
              ref_penalisation_L2_FTD.add(eq);
 
              is >> dim;
              DoubleTab tab = DoubleTab();
              tab.resize(dim);
              for ( int i = 0 ; i<dim ; i++)
                {
                  is >> tab(i);
                  Cerr << "Reading the penalization value along direction " << i+1 <<" :" << tab(i) << finl;
                }
              tab_penalisation_L2_FTD.add(tab);
              is >> mot ;
            }
        }
      else
        {
          Cerr << "Error while reading the parameters of the L2 penalization " << finl;
          Cerr << "We expected: " << accolade_ouverte << finl;
          exit();
        }
      const Domaine_VF& domaine_vf = ref_cast(Domaine_VF, domaine_dis());
      domaine_vf.domaine().creer_tableau_elements(indic_pena_global);
      domaine_vf.creer_tableau_faces(indic_face_pena_global);
    }
  else
    return Convection_Diffusion_Temperature_base::lire_motcle_non_standard(un_mot,is);
  return 1;
}
 
/*! @brief Associates a physical medium to the equation; the medium is cast to Fluide_base.
 *
 * @param un_milieu the physical medium to associate
 * @throws the medium is not of type Fluide_base
 */
void Convection_Diffusion_Temperature::associer_milieu_base(const Milieu_base& un_milieu)
{
  if (sub_type(Fluide_base,un_milieu)) associer_fluide(ref_cast(Fluide_base, un_milieu));
  else Process::exit(que_suis_je() + " : the fluid " + un_milieu.que_suis_je() + " is not of type Fluide_base!");
}
 
void Convection_Diffusion_Temperature::discretiser()
{
  if (!sub_type(Fluide_reel_base,le_fluide.valeur()))
    if (!le_fluide->has_conductivite() || !le_fluide->has_capacite_calorifique() || !le_fluide->has_beta_t())
      {
        Cerr << "You have not defined the following physical properties of the fluid " << finl;
        Cerr << "needed to solve the energy equation: " << que_suis_je() << finl;
        if  (!le_fluide->has_conductivite()) Cerr << "  Thermal conductivity (lambda)"<< finl;
        if  (!le_fluide->has_capacite_calorifique()) Cerr << "  Specific heat capacity (Cp)"<< finl;
        if  (!le_fluide->has_beta_t()) Cerr << "  Thermal expansion coefficient (beta_th)"<< finl;
        exit();
      }
 
  const Discret_Thyd& dis=ref_cast(Discret_Thyd, discretisation());
  Cerr << "Energy equation discretization" << finl;
  dis.temperature(schema_temps(), domaine_dis(), la_temperature);
  champs_compris_.ajoute_champ(la_temperature);
 
  Equation_base::discretiser();
 
  Cerr << "Convection_Diffusion_Temperature::discretiser() ok" << finl;
}
 
int Convection_Diffusion_Temperature::preparer_calcul()
{
  return Equation_base::preparer_calcul();
}
 
inline int string2int(const char* digit, int& result)
{
  result = 0;
 
  //--- Convert each digit char and add into result.
  while (*digit >= '0' && *digit <='9')
    {
      result = (result * 10) + (*digit - '0');
      digit++;
    }
 
  //--- Check that there were no non-digits at end.
  if (*digit != 0)
    {
      return 0;
    }
 
  return 1;
}
 
void Convection_Diffusion_Temperature::get_noms_champs_postraitables(Noms& nom,Option opt) const
{
  Convection_Diffusion_Temperature_base::get_noms_champs_postraitables(nom,opt);
 
  Noms noms_compris = champs_compris_.liste_noms_compris();
  noms_compris.add("gradient_temperature");
  noms_compris.add("h_echange_");
 
  if (opt == DESCRIPTION)
    Cerr << que_suis_je() << " : " << noms_compris << finl;
  else
    nom.add(noms_compris);
}
 
void Convection_Diffusion_Temperature::creer_champ(const Motcle& motlu)
{
  Convection_Diffusion_Temperature_base::creer_champ(motlu);
 
  if (motlu == "gradient_temperature")
    {
      if (!gradient_temperature)
        {
          const Discret_Thyd& dis=ref_cast(Discret_Thyd, discretisation());
          dis.grad_T(domaine_dis(),domaine_Cl_dis(),la_temperature,gradient_temperature);
          champs_compris_.ajoute_champ(gradient_temperature);
        }
    }
 
  Motcle nom_mot(motlu),temp_mot(nom_mot);
  if (nom_mot.debute_par("H_ECHANGE"))
    {
      if (!h_echange)
        {
          const Discret_Thyd& dis=ref_cast(Discret_Thyd, discretisation());
          temp_mot.suffix("H_ECHANGE_");
          int temperature;
          string2int(temp_mot,temperature);
          dis.h_conv(domaine_dis(),domaine_Cl_dis(),la_temperature,h_echange,nom_mot,temperature);
          champs_compris_.ajoute_champ(h_echange);
        }
    }
}
 
bool Convection_Diffusion_Temperature::has_champ(const Motcle& nom, OBS_PTR(Champ_base)& ref_champ) const
{
  if (nom == "gradient_temperature")
    {
      ref_champ = Convection_Diffusion_Temperature::get_champ(nom);
      return true;
    }
 
  if (h_echange)
    if (nom == h_echange->le_nom())
      {
        ref_champ = Convection_Diffusion_Temperature::get_champ(nom);
        return true;
      }
 
  if (Convection_Diffusion_Temperature_base::has_champ(nom))
    return Convection_Diffusion_Temperature_base::has_champ(nom, ref_champ);
 
  return false; /* nothing found */
}
 
bool Convection_Diffusion_Temperature::has_champ(const Motcle& nom) const
{
  if (nom == "gradient_temperature")
    return true;
 
  if (h_echange)
    if (nom == h_echange->le_nom())
      return true;
 
  if (Convection_Diffusion_Temperature_base::has_champ(nom))
    return true;
 
  return false; /* nothing found */
}
 
const Champ_base& Convection_Diffusion_Temperature::get_champ(const Motcle& nom) const
{
  if (nom == "gradient_temperature")
    {
      double temps_init = schema_temps().temps_init();
      Champ_Fonc_base& ch_gt = ref_cast_non_const(Champ_Fonc_base, gradient_temperature.valeur());
      if (((ch_gt.temps() != la_temperature->temps()) || (ch_gt.temps() == temps_init)) && (la_temperature->mon_equation_non_nul()))
        ch_gt.mettre_a_jour(la_temperature->temps());
      return champs_compris_.get_champ(nom);
    }
 
  if (h_echange)
    if (nom == h_echange->le_nom())
      {
        double temps_init = schema_temps().temps_init();
        Champ_Fonc_base& ch_hconv = ref_cast_non_const(Champ_Fonc_base, h_echange.valeur());
        if (((ch_hconv.temps() != la_temperature->temps()) || (ch_hconv.temps() == temps_init)) && (la_temperature->mon_equation_non_nul()))
          ch_hconv.mettre_a_jour(la_temperature->temps());
        return champs_compris_.get_champ(nom);
      }
 
  return Convection_Diffusion_Temperature_base::get_champ(nom);
}
 
/*! @brief Returns the name of the application domain of the equation.
 *
 * Here "Thermique".
 *
 * @return the name of the application domain of the equation
 */
const Motcle& Convection_Diffusion_Temperature::domaine_application() const
{
  static Motcle domaine = "Thermique";
  return domaine;
}
 
DoubleTab& Convection_Diffusion_Temperature::derivee_en_temps_inco(DoubleTab& derivee)
{
  has_time_factor_ = 1;
  if (!is_penalized)
    {
      return Equation_base::derivee_en_temps_inco(derivee);
    }
  else
    {
      const double rhoCp = get_time_factor();
      // Specific code if temperature equation is penalized
      derivee=inconnue().valeurs();
      //   Setting up the ghost domain method
      mise_en_place_domaine_fantome(derivee);
      DoubleTab& inc=inconnue().valeurs();
      inc = derivee;
 
      DoubleTrav secmem(derivee);
 
      if (nombre_d_operateurs()>1)
        {
          derivee_en_temps_conv(secmem, derivee);
          secmem *= rhoCp;
        }
 
      // IBC and boundary layer transport
      DoubleTrav secmem_conv_vr(derivee);
      //advect the variable with the IBC velocity
      transport_ibc(secmem_conv_vr, derivee);
      //subtract this advection
      secmem -= secmem_conv_vr;
 
      DoubleTab derivee_conv;
      derivee_conv.copy(secmem, RESIZE_OPTIONS::COPY_INIT);
 
      les_sources.ajouter(secmem);
 
      filtrage_si_appart_ibc(derivee_conv,secmem);
      secmem.echange_espace_virtuel();
 
      if (calculate_time_derivative())
        {
          // Store dI/dt(n) = M-1 secmem :
          derivee_en_temps().valeurs()=secmem;
          solveur_masse->appliquer(derivee_en_temps().valeurs());
          schema_temps().modifier_second_membre((*this),secmem); // Change secmem for some schemes (eg: Adams_Bashforth)
        }
 
      solveur_masse->appliquer(secmem);
      Equation_base::Gradient_conjugue_diff_impl(secmem, derivee);
 
      //  temperature penalization
      penalisation_L2(derivee);
      return derivee;
    }
}
 
double Convection_Diffusion_Temperature::get_time_factor() const
{
  return domaine_dis().nb_elem() ? milieu().capacite_calorifique().valeurs()(0, 0) * milieu().masse_volumique().valeurs()(0, 0) : 1.0;
}
 
// Adds contributions from operators and sources
void Convection_Diffusion_Temperature::assembler(Matrice_Morse& matrice, const DoubleTab& inco, DoubleTab& resu)
{
  const double rhoCp = get_time_factor();
 
  // Test of the contribuer_a_avec method
  for (int op=0; op<nombre_d_operateurs(); op++)
    operateur(op).l_op_base().tester_contribuer_a_avec(inco, matrice);
 
  // Contributions from operators and sources:
  // [Vol/dt+A]Inco(n+1)=sum(residual)+Vol/dt*Inco(n)
  // Typically: if Op=flux(Inco) then the implicit matrix A contains a contribution -dflux/dInco
  // Example: convective flux operator in VDF:
  // Op=T*u*S and A=-d(T*u*S)/dT=-u*S
  const Discretisation_base::type_calcul_du_residu& type_codage=probleme().discretisation().codage_du_calcul_du_residu();
  if (type_codage==Discretisation_base::VIA_CONTRIBUER_AU_SECOND_MEMBRE)
    {
      if ( probleme().discretisation().que_suis_je() == "EF")
        {
          // Compute sum(residual) via contribuer_au_second_membre (typically for non-implicitized BCs)
          // This approach requires coding 3 methods (contribuer_a_avec, contribuer_au_second_membre, and ajouter for the explicit scheme)
          sources().contribuer_a_avec(inco,matrice);
          statistics().end_count(STD_COUNTERS::matrix_assembly,0,0);
          sources().ajouter(resu);
          statistics().begin_count(STD_COUNTERS::matrix_assembly,statistics().get_last_opened_counter_level()+1);
          matrice.ajouter_multvect(inco, resu); // Add source residual first
          for (int op = 0; op < nombre_d_operateurs(); op++)
            {
              operateur(op).l_op_base().contribuer_a_avec(inco, matrice);
              operateur(op).l_op_base().contribuer_au_second_membre(resu);
            }
        }
      else
        {
          Cerr << "VIA_CONTRIBUER_AU_SECOND_MEMBRE not implemented for " << que_suis_je() << ":assembler" << finl;
          Cerr << "with discretization " <<  probleme().discretisation().que_suis_je() << "" << finl;
          Process::exit();
        }
      // // Compute sum(residual) via contribuer_au_second_membre (typically for non-implicitized BCs)
      // // This approach requires coding 3 methods (contribuer_a_avec, contribuer_au_second_membre, and ajouter for the explicit scheme)
      // sources().contribuer_a_avec(inco,matrice);
      // sources().ajouter(resu);
      // matrice.ajouter_multvect(inco, resu); // Add source residual first
      // for (int op = 0; op < nombre_d_operateurs(); op++)
      //   {
      //     operateur(op).l_op_base().contribuer_a_avec(inco, matrice);
      //     operateur(op).l_op_base().contribuer_au_second_membre(resu);
      //   }
    }
  else if (type_codage==Discretisation_base::VIA_AJOUTER)
    {
      // Compute sum(residual) via sum(operator)+sources+A*Inco(n)
      // This approach requires coding only two methods (contribuer_a_avec and ajouter)
      // Slightly more expensive computationally but less code to write/maintain
      for (int op=0; op<nombre_d_operateurs(); op++)
        {
          Matrice_Morse mat(matrice);
          mat.get_set_coeff() = 0.0;
          operateur(op).l_op_base().contribuer_a_avec(inco, mat);
          if (op == 1) mat *= rhoCp; // the time derivative is multiplied by rhoCp for convection
          matrice += mat;
          statistics().end_count(STD_COUNTERS::matrix_assembly,0,0);
          {
            DoubleTab resu_tmp(resu);
            resu_tmp = 0.;
            operateur(op).ajouter(inco, resu_tmp);
            if (op == 1) resu_tmp *= rhoCp;
            resu += resu_tmp;
          }
          statistics().begin_count(STD_COUNTERS::matrix_assembly,statistics().get_last_opened_counter_level()+1);
        }
 
      sources().contribuer_a_avec(inco,matrice);
      statistics().end_count(STD_COUNTERS::matrix_assembly,0,0);
      sources().ajouter(resu);
      statistics().begin_count(STD_COUNTERS::matrix_assembly,statistics().get_last_opened_counter_level()+1);
      matrice.ajouter_multvect(inco, resu); // Add A*Inco(n)
      // PL (11/04/2018): It would be nice to compute the source contributions first
      // as in the VIA_CONTRIBUER_AU_SECOND_MEMBRE case, but Canal_perio_3D (periodic keps
      // crashes: there is a periodicity error in the source terms of the KEps model...
    }
  else
    {
      Cerr << "Unknown value in Equation_base::assembler for " << (int)type_codage << finl;
      Process::exit();
    }
}
 
void Convection_Diffusion_Temperature::assembler_blocs(matrices_t matrices, DoubleTab& secmem, const tabs_t& semi_impl) const
{
  if (discretisation().is_poly_family() || probleme().que_suis_je().debute_par("Pb_Multiphase"))
    {
      Equation_base::assembler_blocs(matrices, secmem, semi_impl);
      return;
    }
  const double rhoCp = get_time_factor();
  /* mise a zero */
  secmem = 0;
  for (auto &&i_m : matrices) i_m.second->get_set_coeff() = 0;
  matrices_t mats2;
  for (auto &&i_m : matrices)
    {
      Matrice_Morse* mat2 = new Matrice_Morse;
      *mat2=*i_m.second;
      mats2[i_m.first] = mat2;
    }
 
  /* operateurs, sources, masse */
  for (int i = 0; i < nombre_d_operateurs(); i++)
    {
      DoubleTab secmem_tmp(secmem);
      secmem_tmp = 0.;
      for (auto &&i_m : mats2)
        mats2[i_m.first]->get_set_coeff() = 0;
 
      operateur(i).l_op_base().ajouter_blocs(mats2, secmem_tmp, semi_impl);
 
      if (i == 1)
        {
          secmem_tmp *= rhoCp;
          for (auto &&i_m : mats2) *i_m.second *= rhoCp;
        }
 
      for (auto &&i_m : matrices) *i_m.second += *mats2[i_m.first];
      secmem += secmem_tmp;
 
    }
  statistics().end_count(STD_COUNTERS::ajouter_blocs,0,0);
 
  statistics().begin_count(STD_COUNTERS::source_terms,statistics().get_last_opened_counter_level()+1);
  for (int i = 0; i < les_sources.size(); i++)
    les_sources(i)->ajouter_blocs(matrices, secmem, semi_impl);
  statistics().end_count(STD_COUNTERS::source_terms);
 
  statistics().begin_count(STD_COUNTERS::ajouter_blocs,statistics().get_last_opened_counter_level()+1);
 
  const std::string& nom_inco = inconnue().le_nom().getString();
  Matrice_Morse *mat = matrices.count(nom_inco)?matrices.at(nom_inco):nullptr;
  if(mat) mat->ajouter_multvect(inconnue().valeurs(), secmem);
 
  for (auto &&i_m : mats2)
    delete i_m.second;
 
}
int Convection_Diffusion_Temperature::verifier_tag_indicatrice_pena_glob()
{
  //Have the IBCs been modified?
  int maj = 0;
  int tag_all = 0;
  for ( int w = 0; w<ref_penalisation_L2_FTD.size() ; ++w)
    {
      Transport_Interfaces_base& nom_eq = ref_cast(Transport_Interfaces_base,ref_penalisation_L2_FTD[w].valeur());
      tag_all += nom_eq.get_mesh_tag();
    }
  if (tag_all == tag_indic_pena_global) maj = 1;
  return maj;
}
 
int Convection_Diffusion_Temperature::mise_a_jour_tag_indicatrice_pena_glob()
{
  //Have the IBCs been modified?
  int maj = 0;
  int tag_all = 0;
  for ( int w = 0; w<ref_penalisation_L2_FTD.size() ; ++w)
    {
      Transport_Interfaces_base& nom_eq = ref_cast(Transport_Interfaces_base,ref_penalisation_L2_FTD[w].valeur());
      tag_all += nom_eq.get_mesh_tag();
    }
  if (tag_all != tag_indic_pena_global)
    {
      tag_indic_pena_global = tag_all;
      maj = 1;
    }
  return maj;
}
 
void Convection_Diffusion_Temperature::set_indic_pena_globale()
{
  //  Set the global IBC characteristic function (elements and faces)
  int maj = verifier_tag_indicatrice_pena_glob();
  if (maj == 0)
    {
      calcul_indic_pena_global(indic_pena_global, indic_face_pena_global);
      maj = mise_a_jour_tag_indicatrice_pena_glob();
      assert(maj == 1);
    }
}
 
void Convection_Diffusion_Temperature::transport_ibc(DoubleTrav& secmem_conv_vr, DoubleTab& inco_conv_vr)
{
  secmem_conv_vr = 0.;
  const double rhoCp = get_time_factor();
 
  //global IBC characteristic function at the current time
  set_indic_pena_globale();
 
  // compute the thick indicator
  const Domaine_VF& domaine_vf = ref_cast(Domaine_VF, domaine_dis());
  const IntTab& faces_elem = domaine_vf.face_voisins();
  IntTrav indic_pena_global_fat(indic_pena_global);
  for (int i_face = 0; i_face < indic_face_pena_global.size(); i_face++)
    {
      if (indic_face_pena_global(i_face) > 0.)
        {
          const int voisin_0 = faces_elem(i_face, 0);
          if (voisin_0 >= 0) indic_pena_global_fat(voisin_0) = 1;
          const int voisin_1 = faces_elem(i_face, 1);
          if (voisin_1 >= 0) indic_pena_global_fat(voisin_1) = 1;
        }
    }
  // special post-processing for corners
  const IntTab& elem_faces = domaine_vf.elem_faces();
  const int nb_faces_elem = elem_faces.dimension(1);
  const int nb_elem  = indic_pena_global.dimension(0);
  for (int i_elem = 0; i_elem < nb_elem;  i_elem++)
    {
      if (indic_pena_global_fat(i_elem) == 0)
        {
          int coeff = 0;
          for (int i_face = 0; i_face < nb_faces_elem; i_face++)
            {
              const int face = elem_faces(i_elem, i_face);
              const int voisin = faces_elem(face, 0) + faces_elem(face, 1) - i_elem;
              if (voisin >= 0)
                {
                  if (indic_pena_global_fat(voisin) != 0)
                    {
                      ++coeff;
                    }
                }
            }
          if (coeff > 1) indic_pena_global_fat(i_elem) = 1;
        }
    }
 
  // For a local treatment per indicatrice (necessary if using imposed velocity)
  // IntTrav indic_pena_fat(indic_pena_global);
  // for (int w = 0; w < ref_penalisation_L2_FTD.size(); w++)
  //   {
  //     Transport_Interfaces_base & nom_eq = ref_cast(Transport_Interfaces_base,ref_penalisation_L2_FTD[w].valeur());
  //     // on calcul le filtre local a l ibc
  //     const DoubleTab & indicatrice_face = nom_eq.get_compute_indicatrice_faces().valeurs();
  //     indic_pena_fat = 0;
  //     for (int i_face = 0; i_face < indicatrice_face.size(); i_face++)
  //         {
  //           if (indicatrice_face(i_face) > 0.)
  //             {
  //               const int voisin_0 = faces_elem(i_face, 0);
  //               if (voisin_0 >= 0) indic_pena_fat(voisin_0) = 1;
  //               const int voisin_1 = faces_elem(i_face, 1);
  //               if (voisin_1 >= 0) indic_pena_fat(voisin_1) = 1;
  //             }
  //         }
  //     indic_pena_global_fat += indic_pena_fat;
  //   }
 
  // convect the variable field with the fluid velocity for IBC faces
  if ( nombre_d_operateurs() > 1 )
    {
      derivee_en_temps_conv(secmem_conv_vr, inco_conv_vr);
      secmem_conv_vr *= rhoCp;
    }
 
  //keep the contribution only for cells whose fat indicator is != 0
  for (int i_elem = 0; i_elem < nb_elem;  i_elem++)
    {
      if (indic_pena_global_fat(i_elem) == 0) secmem_conv_vr(i_elem) = 0.;
    }
}
 
void Convection_Diffusion_Temperature::mise_en_place_domaine_fantome(DoubleTab& solution)
{
  int maj = verifier_tag_indicatrice_pena_glob();
  if (maj == 1) return;
 
  // safety check
  if (schema_temps().facteur_securite_pas() > 1.0)
    {
      Cerr << "Convection_Diffusion_Temperature::mise_en_place_domaine_fantome safety factor must be <= 1"<<finl;
      Process::exit();
    }
 
  // current global penalized IBC characteristic function
  IntTrav indicatrice_totale;
  IntTrav indicatrice_face_totale ;
  calcul_indic_pena_global(indicatrice_totale,indicatrice_face_totale);
  // current state
  const int nb_elem  = indicatrice_totale.dimension(0);
  int x=0;
  int j=0;
  for (int k_elem =0 ; k_elem <nb_elem; ++k_elem)  ((indicatrice_totale(k_elem)!=0) ? ++x : ++j);
 
  //ghost cell treatment (if any)
  const Domaine_VF& domaine_vf = ref_cast(Domaine_VF, domaine_dis());
  const IntTab& elem_faces = domaine_vf.elem_faces();
  const IntTab& faces_elem = domaine_vf.face_voisins();
  const int nb_faces_elem = elem_faces.dimension(1);
  int k=0; //number of ghost cells ibc -> fluid
  int k_cor=0; //number of ghost cells (ibc -> fluid) corrected
  int k_cor2=0; //number of ghost cells (fluid -> ibc) corrected
  int u=0; //number of fluid -> ibc cells
  for (int i_elem = 0; i_elem < nb_elem; i_elem++)
    {
      //ghost cells (currently fluid, previously ibc) => modify solution
      double somme_inc = 0.;
      double coeff = 0.;
      if (indic_pena_global(i_elem) != 0 && indicatrice_totale(i_elem) == 0)
        {
          ++k;
          // loop over faces to determine the variable value for strictly fluid neighbouring cells
          IntTrav stok_vois(nb_faces_elem);
          for (int i_face = 0; i_face < nb_faces_elem; i_face++)
            {
              const int face = elem_faces(i_elem, i_face);
              const int voisin = faces_elem(face, 0) + faces_elem(face, 1) - i_elem;
              if (voisin >= 0) // the neighbour exists
                {
                  if ((indic_pena_global(voisin) == indicatrice_totale(voisin)) &&  (indic_pena_global(voisin) == indicatrice_totale(i_elem))) // strictly fluid at both time steps
                    {
                      somme_inc += solution(voisin);
                      ++coeff;
                    }
                  else     // save the neighbour index
                    {
                      stok_vois(i_face) = voisin;
                    }
                }
            }
          // If no fluid cells found, check neighbours of neighbours
          if (coeff == 0.)
            {
              for (int i_vois = 0; i_vois < nb_faces_elem; i_vois++)
                {
                  for (int ii_face = 0; ii_face < nb_faces_elem; ii_face++)
                    {
                      const int face_v = elem_faces(stok_vois(i_vois), ii_face);
                      const int voisin_v = faces_elem(face_v, 0) + faces_elem(face_v, 1) - stok_vois(i_vois);
                      if ((voisin_v >= 0) && (voisin_v != i_elem) )// the neighbour exists and is not i_elem
                        {
                          if ((indic_pena_global(voisin_v) == indicatrice_totale(voisin_v)) &&  (indic_pena_global(voisin_v) == indicatrice_totale(i_elem))) // strictly fluid at both time steps
                            {
                              somme_inc += solution(voisin_v);
                              ++coeff;
                            }
                        }
                    }
                }
            }
          // define a new value of the variable for the ghost cell
          if (coeff > 0.)
            {
              solution(i_elem) = somme_inc / coeff;
              k_cor += 1;
            }
        }
      //cells currently ibc and previously fluid
      if (indic_pena_global(i_elem) == 0 && indicatrice_totale(i_elem) != 0)
        {
          ++u;
          // loop over faces to determine the variable value for strictly ibc neighbouring cells
          IntTrav stok_vois(nb_faces_elem);
          for (int i_face = 0; i_face < nb_faces_elem; i_face++)
            {
              const int face = elem_faces(i_elem, i_face);
              const int voisin = faces_elem(face, 0) + faces_elem(face, 1) - i_elem;
              if (voisin >= 0) // the neighbour exists
                {
                  // strictly ibc at both time steps (the same) because solution was modified above
                  if ((indic_pena_global(voisin) == indicatrice_totale(voisin)) && (indic_pena_global(voisin) == indicatrice_totale(i_elem)))
                    {
                      somme_inc += solution(voisin);
                      ++coeff;
                    }
                  else     // save the neighbour index
                    {
                      stok_vois(i_face) = voisin;
                    }
                }
            }
          // If no ibc cells found, check neighbours of neighbours
          if (coeff == 0.)
            {
              for (int i_vois = 0; i_vois < nb_faces_elem; i_vois++)
                {
                  for (int ii_face = 0; ii_face < nb_faces_elem; ii_face++)
                    {
                      const int face_v = elem_faces(stok_vois(i_vois), ii_face);
                      const int voisin_v = faces_elem(face_v, 0) + faces_elem(face_v, 1) - stok_vois(i_vois);
                      if ((voisin_v >= 0) && (voisin_v != i_elem) )// the neighbour exists and is not i_elem
                        {
                          if ((indic_pena_global(voisin_v) == indicatrice_totale(voisin_v)) && (indic_pena_global(voisin_v) == indicatrice_totale(i_elem))) // strictly ibc at both time steps (the same) because solution was modified above
                            {
                              somme_inc += solution(voisin_v);
                              ++coeff;
                            }
                        }
                    }
                }
            }
          // define a new value of the variable for the ghost cell
          if (coeff > 0.)
            {
              solution(i_elem) = somme_inc / coeff;
              k_cor2 += 1;
            }
          else
            {
              solution(i_elem) = tab_penalisation_L2_FTD[indicatrice_totale(i_elem)-1](0);
              k_cor2 += 1;
            }
        }
    }
  solution.echange_espace_virtuel();
  //     Debog::verifier("Convection_Diffusion_Temperature::mise_en_place_domaine_fantome solution ",solution);
  Cout <<"Mise_en_place_domaine_fantome : fluid -> ibc : "<<u<< " elem. (of which "<<k_cor2<<" corrected) and ibc -> fluid : ";
  Cout<<k<<" elem.(of which "<<k_cor<<" corrected)"<<finl;
  Cout <<"             nb ibc elem. : "<<x<<" and nb fluid elem. : "<<j<<" out of a total of : "<<nb_elem<<finl;
  if ((u != k_cor2) || (k !=k_cor))
    {
      Cerr <<"Mise_en_place_domaine_fantome : Not all ghost cells have been corrected"<<finl;
      exit();
    }
 
  indic_pena_global = indicatrice_totale;
  indic_face_pena_global = indicatrice_face_totale;
  maj = mise_a_jour_tag_indicatrice_pena_glob();
  assert(maj == 1);
}
 
DoubleTab& Convection_Diffusion_Temperature::filtrage_si_appart_ibc(DoubleTab& u_conv, DoubleTab& u)
{
  // compute the global IBC characteristic function for the neighbourhood temperature T_voisinage
  set_indic_pena_globale();
 
  for (int j = 0 ; j< indic_pena_global.size() ; ++j)
    {
      if ( indic_pena_global(j) != 0) u(j) = u_conv(j);
    }
  u.echange_espace_virtuel();
  //   Debog::verifier("Convection_Diffusion_Temperature::filtrage_si_appart_ibc u ",u);
  return u;
}
 
void Convection_Diffusion_Temperature::calcul_indic_pena_global(IntTab& indicatrice_totale, IntTab& indicatrice_face_totale)
{
  // global penalized IBC characteristic function
  indicatrice_totale.copy(indic_pena_global);
  indicatrice_face_totale.copy(indic_face_pena_global);
  indicatrice_totale = 0;
  indicatrice_face_totale = 0;
  int x=0;
  int xf=0;
  const int nb_elem  = indicatrice_totale.dimension_tot(0);
  const int nfaces = indicatrice_face_totale.dimension_tot(0);
  // loop over ibcs
  for ( int w = 0 ; w<ref_penalisation_L2_FTD.size() ; ++w)
    {
      Transport_Interfaces_base& nom_eq = ref_cast(Transport_Interfaces_base,ref_penalisation_L2_FTD[w].valeur());
      nom_eq.check_indicatrice_is_up_to_date();
      const DoubleTab& indicatrice = nom_eq.get_indicatrice().valeurs();
      // characteristic function (0 or 1) for all IBCs
      for (int k_elem =0 ; k_elem <nb_elem; ++k_elem)
        {
          if (indicatrice(k_elem)>0.)
            {
              if(indicatrice_totale(k_elem) != 0)
                {
                  Cerr<<"calcul_indic_pena_global : Warning: the elements of the indicator functions overlap "<<finl;
                  exit();
                }
              indicatrice_totale(k_elem) = w+1;
              ++x;
            }
        }
      // characteristic function (IBC number) for all IBCs
      const Domaine_VF& domaine_vf = ref_cast(Domaine_VF, domaine_dis());
      const IntTab& face_voisins = domaine_vf.face_voisins();
      for (int i = 0; i < nfaces ; i++)
        {
          const int elem0 = face_voisins(i, 0);
          const int elem1 = face_voisins(i, 1);
          double indi = 0.;
          if (elem0 >= 0) indi = indicatrice(elem0);
          if (elem1 >= 0) indi += indicatrice(elem1);
          if(indi > 0.)
            {
              if(indicatrice_face_totale(i) !=  0)
                {
                  Cerr<<"calcul_indic_pena_global : Warning: the faces of the indicator functions overlap "<<finl;
                  exit();
                }
              indicatrice_face_totale(i) = w+1;
              ++xf;
            }
        }
    }
  Cout<<"Calcul_indic_pena_global : Number of elements of the global IBC indicator: "<<x;
  Cout<<"  and number of faces: "<<xf<<finl;
  indicatrice_totale.echange_espace_virtuel();
  indicatrice_face_totale.echange_espace_virtuel();
  //      Debog::verifier("Convection_Diffusion_Temperature::calcul_indic_pena_global indicatrice_totale ",indicatrice_totale);
  //      Debog::verifier("Convection_Diffusion_Temperature::calcul_indic_pena_global indicatrice_face_totale ",indicatrice_face_totale);
}
 
DoubleTab& Convection_Diffusion_Temperature::penalisation_L2(DoubleTab& u)
{
  // garde-fou
  assert(ref_penalisation_L2_FTD.size() != 0);
  assert(ref_penalisation_L2_FTD.size() == tab_penalisation_L2_FTD.size());
 
  const double dt = schema_temps().pas_de_temps();
  DoubleTab u_old(u);
 
  // compute the global IBC characteristic function for the neighbourhood temperature T_voisinage
  set_indic_pena_globale();
 
  //compute T_voisinage for all elements
  const Domaine_VF& domaine_vf = ref_cast(Domaine_VF, domaine_dis());
  const IntTab& elem_faces = domaine_vf.elem_faces();
  const IntTab& faces_elem = domaine_vf.face_voisins();
  const int nb_faces_elem = elem_faces.dimension(1);
  const DoubleTab& inc=inconnue().valeurs();
  // the unknown must be scalar
  assert(inc.line_size() == 1);
  DoubleTrav t_voisinage(inc);
  DoubleTrav u_voisinage(u);
  const int nb_elem  = u.dimension_tot(0);
  for (int i_elem = 0; i_elem < nb_elem; i_elem++)
    {
      double somme = 0.;
      double somme_deriv = 0.;
      double coeff = 0.;
 
      if( indic_pena_global(i_elem) != 0)
        {
          // loop over faces of the element
          for (int i_face = 0; i_face < nb_faces_elem; i_face++)
            {
              const int face = elem_faces(i_elem, i_face);
              const int voisin = faces_elem(face, 0) + faces_elem(face, 1) - i_elem;
              if (voisin >= 0)
                {
                  if (indic_pena_global(voisin) == 0)
                    {
                      somme += inc(voisin);
                      somme_deriv += u(voisin);
                      coeff += 1.;
                    }
                }
            }
        }
 
      if (coeff > 0.)
        {
          t_voisinage(i_elem) = somme / coeff;
          u_voisinage(i_elem) = somme_deriv / coeff;
        }
      else
        {
          // case: no fluid neighbour
          t_voisinage(i_elem) = inc(i_elem);
          u_voisinage(i_elem) = u(i_elem);
        }
    }
  t_voisinage.echange_espace_virtuel();
  u_voisinage.echange_espace_virtuel();
 
  //determination of the global L2 penalization field
  DoubleTrav pena_glob(inc); //initialise a zero
  DoubleTrav tab_(inc); //initialise a zero
  DoubleTrav denom(inc); //initialise a zero
 
  //loop over IBCs
  for ( int i = 0 ; i < ref_penalisation_L2_FTD.size() ; ++i)
    {
      Transport_Interfaces_base& nom_eq = ref_cast(Transport_Interfaces_base,ref_penalisation_L2_FTD[i].valeur());
      const DoubleTab& tab = tab_penalisation_L2_FTD[i];
      //check dimensions: tab t obstacle = nb dimensions of the unknown
      if ( tab.dimension(0) != inc.line_size() )
        {
          Cerr << " penalisation_L2: The imposed and computed temperature fields do not have the same dimensions" <<finl;
          Cerr << " Dimension of the imposed temperature field "<< tab.dimension(0) << finl;
          Cerr << " Dimension of the temperature field computed by TRUST "<<inc.line_size() <<finl;
          Cerr << " Number of elements stored in the unknown " <<inc.size() <<finl;
          Process::exit();
        }
      nom_eq.check_indicatrice_is_up_to_date();
      const DoubleTab& indicatrice = nom_eq.get_indicatrice().valeurs();
 
      DoubleTrav pena_loc(indicatrice);
      // Determine the penalization term for the current IBC
      // loop over elements
      for (int j = 0; j<nb_elem ; ++j)
        {
          //compute the characteristic function Ksi for the current IBC -> pena_loc
          if (indicatrice(j) > 0.)
            {
              double coeff_p=1.;
              pena_loc(j) = coeff_p;
              //compute the denominator of the penalization term ( Sigma_ibc Ksi_ibc(j) )
              denom(j) += coeff_p;
 
              //determine the value to impose T_ref
              double tref_j = 0.;
              switch (choix_pena)
                {
                case 0 :
                  // approximation of T_ref by T_ref = T_imp
                  tref_j = tab(0);
                  break;
                case 1 :
                  // approximation of T_ref by T_ref = T+dt*derivative
                  tref_j = inc(j) + dt*u_old(j);
                  break;
                case 2 :
                  // approximation of T_imp by T_ref= tau*T_imp + (1-tau)*T_neighbour
                  // approximation of T_neighbour by T_neighbour= 1/n * sum(1,n)(T_surrounding+dt*derivative_surrounding)
                  tref_j = (indicatrice(j)*tab(0)+(1.-indicatrice(j))*(t_voisinage(j)+dt*u_voisinage(j)));
                  break;
                case 3 :
                  // approximation of T_imp by T_ref= [(1-tau)*T_neighbour + (1/2)*T_imp] / (3/2 -tau)
                  // approximation of T_neighbour by T_neighbour= 1/n * sum(1,n)(T_surrounding+dt*derivative_surrounding)
                  tref_j = ((1.-indicatrice(j))*(t_voisinage(j)+dt*u_voisinage(j))+0.5*tab(0))/(1.5-indicatrice(j));
                  break;
                case 4 :
                  // approximation of T_imp by T_ref= T_imp if tau=1, else T_ref=T_neighbour
                  if (indicatrice(j) == 1.)
                    {
                      tref_j = tab(0);
                    }
                  else
                    {
                      tref_j = (t_voisinage(j)+dt*u_voisinage(j));
                    }
                  break;
                default:
                  Cerr << "Penalisation_L2: invalid choix_pena "<< choix_pena <<finl;
                  exit();
                  break;
                }
              // store Sigma_ibc Ksi_ibc(j) Tref_ibc for post-processing
              tab_(j) += tref_j;
              // penalisation_global = Sigma_ibc Ksi_ibc(j)*(Tref_ibc - inc) /  eta;
              pena_loc(j) *= (tref_j-inc(j));
              pena_loc(j) /= eta;
              pena_glob(j)  += pena_loc(j);
            }
        }
    }
 
  //  Apply the L2 penalization of the time derivative: (dT_t_* + pena_glob) / ( 1 + (dt*Sigma_ibc Ksi_ibc(j)) / eta )
  // Tref is the arithmetic mean of Tref_j when 0<eta<<1
  for (int j = 0 ; j<nb_elem ; ++j)
    {
      u(j) = u_old(j) + pena_glob(j);
      u(j) /= (1. + denom(j)*dt/eta);
    }
  u.echange_espace_virtuel();
 
  //  Save fluxes
  ecrire_fichier_pena_th(u_old,u,tab_,denom);
 
  //   Debog::verifier("Convection_Diffusion_Temperature::penalisation_L2 u ",u);
  return u;
}
 
void ouvrir_fichier_pena_th(SFichier& os, const Nom& type, const int flag, const Transport_Interfaces_base& equation)
{
  // flag nul on n'ouvre pas le fichier
  if (flag==0)
    return ;
  Nom fichier=Objet_U::nom_du_cas();
  fichier+="_";
  fichier+=type;
  fichier+="_";
  fichier+=equation.le_nom();
  fichier+=".out";
  const Schema_Temps_base& sch=equation.probleme().schema_temps();
  const int precision=sch.precision_impr();
  // On cree le fichier a la premiere impression avec l'en tete
  if (sch.nb_impr()==1 && !equation.probleme().reprise_effectuee())
    {
      os.ouvrir(fichier);
      SFichier& fic=os;
      Nom espace="\t\t";
      fic << (Nom)"# Impression puissance thermique";
      fic << " interface " << equation.le_nom();
      fic << " en W" << "." << finl;
      fic << finl << "# Temps";
 
      Nom ch=espace;
      ch+="P";
      fic << ch << finl;
    }
  // Sinon on l'ouvre
  else
    {
      os.ouvrir(fichier,ios::app);
    }
  os.precision(precision);
  os.setf(ios::scientific);
}
 
void Convection_Diffusion_Temperature::ecrire_fichier_pena_th(DoubleTab& u_old, DoubleTab& u, DoubleTab& tref, DoubleTab& denom)
{
 
  if (le_schema_en_temps->limpr())
    {
 
      const Domaine_VF& domaine_vf = ref_cast(Domaine_VF, domaine_dis());
      const DoubleVect& vol_maille = domaine_vf.volumes();
      const Fluide_base& fluide_inc = ref_cast(Fluide_base, milieu());
      const DoubleTab& tab_rho = fluide_inc.masse_volumique().valeurs();
      const double rho = tab_rho(0,0);
      const DoubleTab& tab_cp = fluide_inc.capacite_calorifique().valeurs();
      const double cp = tab_cp(0,0);
 
      //  Method to compute the total flux on the IBCs
 
      //Compute T_n+1 after penalization
      const DoubleTab& inc=inconnue().valeurs();
      const int nb_elem  = u.dimension(0);
      const double dt = schema_temps().pas_de_temps();
      DoubleTab tkp1(inc);
      for (int k = 0 ; k<nb_elem ; ++k) tkp1(k) += u(k)*dt;
 
      //Terme de penalisation:       Sigma_ibc Ksi_ibc * ( T_n+1 - Tref_ibc) / eta
      DoubleTrav F(inc);
      for (int j = 0; j<nb_elem; ++j)  F(j) = (tref(j)-denom(j)*tkp1(j))/eta;
 
      // loop over IBCs
 
      DoubleTrav filtre(nb_elem);
      DoubleTrav filtre_glob(nb_elem);
      const double temps_flux = le_schema_en_temps->temps_courant();
      Nom espace=" \t";
 
      for ( int i = 0 ; i < ref_penalisation_L2_FTD.size() ; ++i) //loop over the number of IBCs
        {
          Transport_Interfaces_base& nom_eq = ref_cast(Transport_Interfaces_base,ref_penalisation_L2_FTD[i].valeur());
          const DoubleTab& indicatrice = nom_eq.get_indicatrice().valeurs();
          double Flux_pena= 0.;
          double Flux_pena_old= 0.;
          double Flux_temp_interne= 0.;
          //Filter: ibc (1) / non-ibc (0)
          filtre = 0.;
          for ( int j = 0 ; j<nb_elem ; ++j)
            {
              if (indicatrice(j) > 0.)
                {
                  filtre(j) = 1.;
                  filtre_glob(j)=1.;
                }
              //bilans
              //Note: rho and cp are constant
              Flux_pena_old += u_old(j) * filtre(j) *rho *cp *vol_maille(j);
              Flux_pena += F(j) * filtre(j) * vol_maille(j) *rho *cp;
              Flux_temp_interne += filtre(j)*rho*cp*vol_maille(j)*u(j);
            }
 
          // Sum over all parallel processes
          Flux_pena_old = mp_sum(Flux_pena_old);
          Flux_pena = mp_sum(Flux_pena);
          Flux_temp_interne = mp_sum(Flux_temp_interne);
 
          //write output
          if (Process::je_suis_maitre())
            {
              SFichier FTE;
              ouvrir_fichier_pena_th(FTE,"Puissance_penalisation_thermique_ibc",1,nom_eq);
              FTE << temps_flux;
              FTE << espace << Flux_pena;
              FTE << finl;
 
              SFichier FTI;
              ouvrir_fichier_pena_th(FTI,"Puissance_thermique_ibc",1,nom_eq);
              FTI << temps_flux;
              FTI << espace << Flux_pena_old;
              FTI << finl;
 
              SFichier FTTI;
              ouvrir_fichier_pena_th(FTTI,"Derivee_temps_temperature_ibc",1,nom_eq);
              FTTI << temps_flux;
              FTTI << espace << Flux_temp_interne;
              FTTI << finl;
            }
        }
 
      //budget for Derivee_temps_temperature_fluide
      //Note: rho and cp are constant
      double Flux_temp_externe= 0.;
      for (int j = 0; j<nb_elem; ++j)
        {
          Flux_temp_externe += (1.-filtre_glob(j))  *rho *cp *vol_maille(j) * u(j); //fluid
        }
 
      // Sum contributions from all processors in parallel
      Flux_temp_externe = mp_sum(Flux_temp_externe); //fluid
 
      //write output
      if (Process::je_suis_maitre())
        {
          SFichier fichier;
          Nom fichier_nom=Objet_U::nom_du_cas();
          fichier_nom+="_";
          fichier_nom+="Derivee_temps_temperature_fluide.out";
          if (le_schema_en_temps->nb_impr()==1 && !probleme().reprise_effectuee())
            {
              fichier.ouvrir(fichier_nom);
            }
          else
            {
              fichier.ouvrir(fichier_nom,ios::app);
            }
          //if (fichier)
          {
            (fichier) <<temps_flux<<" \t"<<Flux_temp_externe<<finl;
 
          }
        }
    }
}
 
void Convection_Diffusion_Temperature::calculer_rho_cp_T(const Objet_U& obj, DoubleTab& val, DoubleTab& bval, tabs_t& deriv)
{
  const Equation_base& eqn = ref_cast(Equation_base, obj);
  const Fluide_base& fl = ref_cast(Fluide_base, eqn.milieu());
  const Champ_Inc_base& ch_T = eqn.inconnue();
  const Champ_base& ch_rho = fl.masse_volumique();
  assert(sub_type(Champ_Uniforme, ch_rho));
  const Champ_Don_base& ch_cp = fl.capacite_calorifique();
  const DoubleTab& cp = fl.capacite_calorifique().valeurs(), &rho = ch_rho.valeurs(), &T = ch_T.valeurs();
 
  /* valeurs du champ */
  const int N = val.line_size(), Nl = val.dimension_tot(0), cCp = sub_type(Champ_Uniforme, ch_cp);
  for (int i = 0; i < Nl; i++)
    for (int n = 0; n < N; n++) val(i, n) = rho(0, n) * cp(!cCp * i, n) * T(i, n);
 
  /* valeur_aux_bords can only be used if ch_rho has a domaine_dis_base */
  DoubleTab b_cp = cCp ? cp : ch_cp.valeur_aux_bords(), b_T = ch_T.valeur_aux_bords();
  int Nb = b_T.dimension_tot(0);
  // rho is assumed to be a champ_uniforme: use the field array directly
  for (int i = 0; i < Nb; i++)
    for (int n = 0; n < N; n++) bval(i, n) = b_cp(!cCp * i, n) * rho(0, n) * b_T(i, n);
 
  DoubleTab& d_T = deriv["temperature"];
  d_T.resize(Nl, N);
  for (int i = 0; i < Nl; i++)
    for (int n = 0; n < N; n++) d_T(i, n) = rho(0, n) * cp(!cCp * i, n);
 
}