Source_PDF_base.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 <Domaine_Cl_dis_base.h>
#include <Schema_Temps_base.h>
#include <Source_PDF_base.h>
#include <Equation_base.h>
#include <Probleme_base.h>
#include <Type_info.h>
#include <SFichier.h>
#include <Param.h>
#include <Interpolation_IBM_elem_fluid.h>
#include <Interpolation_IBM_mean_gradient.h>
#include <Interpolation_IBM_hybrid.h>
#include <Interpolation_IBM_power_law_tbl.h>
#include <Interpolation_IBM_power_law_tbl_u_star.h>
#include <Interpolation_IBM_thermal_wall_law.h>
 
Implemente_base(Source_PDF_base,"Source_PDF_base",Source_dep_inco_base);
// XD source_pdf source_base source_pdf BRACE Source term for Penalised Direct Forcing (PDF) method.
 
Entree& Source_PDF_base::readOn(Entree& s)
{
  Param param(que_suis_je());
  param.ajouter("prepro_ibm", &prepro_lu_,Param::OPTIONAL); // XD_ADD_P Prepro_IBM_base
  // XD_CONT to realise the PDF IBM preprocessing
  param.ajouter("aire", &champ_aire_lu_,Param::OPTIONAL); // XD_ADD_P field_base
  // XD_CONT volumic field: a boolean for the cell (0 or 1) indicating if the obstacle is in the cell
  param.ajouter_flag("get_aire_from_prepro", &aire_from_prepro_); // XD_ADD_P rien
  // XD_CONT get aire IBM from prepro.
  param.ajouter("barycentre", &champ_barycentre_lu_,Param::OPTIONAL); // XD_ADD_P field_base
  // XD_CONT volumic field with 3 components representing the face barycenters
  param.ajouter_flag("get_barycenter_from_prepro", &barycentre_from_prepro_); // XD_ADD_P rien
  // XD_CONT get aire IBM from prepro.
  param.ajouter("rotation", &champ_rotation_lu_,Param::OPTIONAL); // XD_ADD_P field_base
  // XD_CONT volumic field with 9 components representing the change of basis on cells (local to global). Used for
  // XD_CONT rotating cases for example.
  param.ajouter_flag("get_rotation_from_prepro", &rotation_from_prepro_); // XD_ADD_P rien
  // XD_CONT get aire IBM from prepro.
  param.ajouter_flag("transpose_rotation", &transpose_rotation_); // XD_ADD_P rien
  // XD_CONT whether to transpose the basis change matrix.
  param.ajouter("modele",&modele_lu_,Param::REQUIRED);   // XD_ADD_P bloc_pdf_model
  // XD_CONT model used for the Penalized Direct Forcing
  temps_relax_ = modele_lu_.temps_relax_;
  echelle_relax_ =  modele_lu_.echelle_relax_;
  param.ajouter("interpolation",&interpolation_lue_,Param::OPTIONAL); // XD_ADD_P interpolation_ibm_base
  // XD_CONT interpolation method
 
  param.lire_avec_accolades(s);
  if (interpolation_lue_)
    {
      if (!(interpolation_lue_->que_suis_je() == "Interpolation_IBM_aucune"))
        {
          interpolation_bool_ = true;
          if ((interpolation_lue_->que_suis_je() == "Interpolation_IBM_power_law_tbl") || (interpolation_lue_->que_suis_je() == "Interpolation_IBM_power_law_tbl_u_star") ) imm_wall_law_ = true;
          interpolation_lue_->set_source(*this);
        }
    }
 
  if (!equation().probleme().que_suis_je().contient("_IBM"))
    {
      Cerr << "Error in the equation " << equation().que_suis_je() << " !" << finl;
      Cerr << "You are using a " << que_suis_je() << " source term but your problem is not an IBM problem !" << finl;
      Cerr << "Fix your data set and use one of the following problems : " << finl;
 
      Noms noms_pb;
      Type_info::les_sous_types(Nom("Probleme_base"), noms_pb);
 
      for (auto &itr : noms_pb)
        if (itr.contient("_IBM"))
          Cerr << "  - " << itr << finl;
 
      Cerr << "  - Pb_Couple_Optimisation_IBM" << finl;
 
      Process::exit();
    }
 
  return s;
}
 
Sortie& Source_PDF_base::printOn(Sortie& s ) const
{
  return s << que_suis_je() ;
}
 
void Source_PDF_base::associer_pb(const Probleme_base& pb)
{
  if (prepro_lu_)
    {
      prepro_lu_->associer_pb(pb);
 
      // Optional verification
      if( prepro_lu_->verify_results_prepro() == 1) verify_results_prepro();
    }
 
  const Domaine_dis_base& le_dom_dis = pb.domaine_dis();
 
  // Rotation Matrix
  int dim_esp = Objet_U::dimension;
  assert(dim_esp==3);
  int nb_comp=dim_esp*dim_esp;
  Noms nom_c(nb_comp);
  Noms unites(nb_comp);
  pb.discretisation().discretiser_champ("champ_elem",le_dom_dis,vectoriel,nom_c,unites,nb_comp,0.,champ_rotation_);
  if (prepro_lu_ && rotation_from_prepro_)
    {
      champ_rotation_->valeurs() = prepro_lu_->get_champ_rotation();
    }
  else
    {
      if (champ_rotation_lu_)
        {
          champ_rotation_->affecter(champ_rotation_lu_);
        }
      else
        {
          Cerr<<"Source_PDF_base::associer_pb: rotation term is missing "<<finl;
          exit();
        }
    }
 
  // Barycenter
  int nb_comp0=dim_esp;
  Noms nom_c0(nb_comp0);
  Noms unites0(nb_comp0);
  pb.discretisation().discretiser_champ("champ_elem",le_dom_dis,vectoriel,nom_c0,unites0,nb_comp0,0.,champ_barycentre_);
  if (prepro_lu_ && barycentre_from_prepro_)
    {
      champ_barycentre_->valeurs() = prepro_lu_->get_champ_barycentre();
    }
  else
    {
      if (champ_barycentre_lu_)
        {
          champ_barycentre_->affecter(champ_barycentre_lu_);
        }
      // else
      // barycenters not required without prepro_ibm
      //   {
      //     Cerr<<"Source_PDF_base::associer_pb: barycenter term is missing "<<finl;
      //     exit();
      //   }
    }
 
  // Field for PDF source term
  if (equation().discretisation().is_ef())
    pb.discretisation().discretiser_champ("champ_sommets",le_dom_dis,"","",1,0., champ_nodal_);
  else
    pb.discretisation().discretiser_champ("champ_face",le_dom_dis,"","",1,0., champ_nodal_);
 
  // Area
  pb.discretisation().discretiser_champ("champ_elem",le_dom_dis ,"aire","m-1",1,0., champ_aire_);
  if (prepro_lu_ && aire_from_prepro_)
    {
      champ_aire_->valeurs() = prepro_lu_->get_champ_aire();
    }
  else
    {
      if (champ_aire_lu_)
        {
          champ_aire_->affecter(champ_aire_lu_);
        }
      else
        {
          Cerr<<"Source_PDF_base::associer_pb: aire term is missing "<<finl;
          exit();
        }
    }
 
  // Imposed Shape IB Velocity
  if (get_modele().get_PDF_mobile())
    {
      nb_comp=dim_esp;
      Noms nom_c11(nb_comp);
      Noms unites11(nb_comp);
      pb.discretisation().discretiser_champ("champ_elem",le_dom_dis,vectoriel,nom_c11,unites11,nb_comp,0.,modele_lu_.vitesse_shape_IBM_);
    }
 
  // Variable Imposed on IB
  const DoubleTab& variable=equation().inconnue().valeurs();
  Motcle directive("temperature");
  if (equation().inconnue().is_vectorial()) directive="vitesse";
  nb_comp=variable.dimension(1);
  assert(nb_comp==modele_lu_.dim_variable_);
  Nom nom_c1(equation().inconnue().le_nom());
  Nom unites1(equation().inconnue().unites()[0]);
  type_variable_imposee_ = modele_lu_.type_variable_imposee_;
  pb.discretisation().discretiser_champ(directive,le_dom_dis ,nom_c1,unites1,nb_comp,0.,modele_lu_.variable_imposee_);
 
  if (interpolation_bool_)
    {
      Domaine_dis_base& le_dom_dis_base = ref_cast_non_const(Domaine_dis_base,le_dom_dis);
      interpolation_lue_->discretise(pb.discretisation(),le_dom_dis_base);
    }
 
  // set variable imposee
  set_variable_imposee();
 
  // Rho
  pb.discretisation().discretiser_champ("champ_elem",le_dom_dis ,"rho","kg.m-3",1,0., champ_rho_);
  Source_PDF_base::updateChampRho();
 
  // Relaxation
  temps_relax_ = modele_lu_.temps_relax_;
  echelle_relax_ = modele_lu_.echelle_relax_;
 
  matrice_pression_variable_bool_ = false;
  if (temps_relax_ != 1.0e+12) matrice_pression_variable_bool_ = true;
}
 
void Source_PDF_base::set_variable_imposee()
{
  int dim_esp = Objet_U::dimension;
  const Probleme_base& pb = equation().probleme();
  const Domaine_dis_base& le_dom_dis = pb.domaine_dis();
 
  // Transpose Rotation Matrix
  if (transpose_rotation_)
    {
      DoubleTab& val=champ_rotation_->valeurs();
 
      int nb_case=val.dimension_tot(0);
      for (int ele=0; ele<nb_case; ele++)
        for (int k=0; k<dim_esp; k++)
          for (int i=k+1; i<dim_esp; i++)
            {
              double tmp=val(ele,3*k+i);
              val(ele,3*k+i)=val(ele,3*i+k);
              val(ele,3*i+k)=tmp;
            }
    }
 
  // Variable Imposed on IB
  const DoubleTab& variable=equation().inconnue().valeurs();
  int nb_comp=variable.dimension(1);
 
  if (interpolation_bool_) // interpolation
    {
      if (type_variable_imposee_ == 1) // fonction xyzt
        {
          if (interpolation_lue_->que_suis_je() == "Interpolation_IBM_gradient_moyen" || interpolation_lue_->que_suis_je() == "Interpolation_IBM_power_law_tbl_u_star")
            {
              const Interpolation_IBM_mean_gradient& interp = ref_cast(Interpolation_IBM_mean_gradient,interpolation_lue_.valeur());
              this->compute_variable_imposee_projete(interp.solid_elems_->valeurs(), interp.solid_points_->valeurs(), -2.0, 1e-6);
            }
          else
            {
              const Interpolation_IBM_elem_fluid& interp = ref_cast(Interpolation_IBM_elem_fluid,interpolation_lue_.valeur());
              compute_variable_imposee_projete(interp.fluid_elems_->valeurs(), interp.solid_points_->valeurs(), -2.0, 1e-6);
            }
        }
      else // data
        {
          modele_lu_.variable_imposee_->affecter(modele_lu_.variable_imposee_lu_);
        }
    }
  else // no interpolation
    {
      if (type_variable_imposee_ != 1) // data
        {
          modele_lu_.variable_imposee_->affecter(modele_lu_.variable_imposee_lu_);
        }
      else // fonction xyzt
        {
          const DoubleTab& coords = le_dom_dis.domaine().coord_sommets();
          Domaine_VF& le_dom_VF = ref_cast_non_const(Domaine_VF, pb.domaine_dis());
          modele_lu_.affecter_variable_imposee(le_dom_VF, coords);
        }
    }
  if (modele_lu_.local_ == 1)
    {
      if (nb_comp != dim_esp)
        {
          Cerr<<"Source_PDF_base::calcul_variable_imposee: variable dimension different from "<<dim_esp<<"! "<<finl;
          abort();
        }
      rotate_imposed_velocity(modele_lu_.variable_imposee_->valeurs());
    }
 
  variable_imposee_ = modele_lu_.variable_imposee_->valeurs();
 
  compute_indicateur_nodal_champ_aire();
}
 
void Source_PDF_base::set_fields_to_prepro(Prepro_IBM_base& un_prepro)
{
  un_prepro.set_champ_aire(champ_aire_->valeurs());
  un_prepro.set_champ_rotation(champ_rotation_->valeurs());
  un_prepro.set_champ_barycentre(champ_barycentre_->valeurs());
}
 
void Source_PDF_base::get_fields_from_prepro(Prepro_IBM_base& un_prepro)
{
  champ_aire_->valeurs() = un_prepro.get_champ_aire();
  champ_rotation_->valeurs() = un_prepro.get_champ_rotation();
  champ_barycentre_->valeurs() = un_prepro.get_champ_barycentre();
}
 
void Source_PDF_base::compute_NeighNode_IBM_elem(DoubleTab& aireArray, IntLists& elem_voisins, bool all_elem_vois)
{
  const Probleme_base& pb = equation().probleme();
  const Domaine_dis_base& le_dom_dis = pb.domaine_dis();
  const Domaine_VF& the_dom_VF = ref_cast(Domaine_VF,le_dom_dis);
  const Domaine& le_dom =  le_dom_dis.domaine();
  int nb_elem_tot = le_dom.nb_elem_tot();
  int nb_som_elem = le_dom.nb_som_elem();
  int nb_faces_elem = le_dom.nb_faces_elem();
  const IntTab& face_voisins = the_dom_VF.face_voisins();
  const IntTab& elems = le_dom.les_elems() ;
  const IntTab& elem_face = the_dom_VF.elem_faces();
 
  // Clear the Lists
  assert (elem_voisins.size() == nb_elem_tot);
  for (int num_elem = 0; num_elem < nb_elem_tot; num_elem++) elem_voisins[num_elem].vide();
 
  for (int num_elem = 0; num_elem < nb_elem_tot; num_elem++) // includes ghost elements
    {
      if (aireArray(num_elem) > 0.)
        {
          if (!elem_voisins[num_elem].est_vide()) // Liste non vide !
            {
              Cerr<<"IBM : Source_PDF_base::compute_NeighNode_IBM_elem:  Neighbour liste is not empty for element "<<num_elem<<" : Abort"<<finl;
              exit();
            }
          // Cerr << "******Elem = "<<num_elem<<finl;
 
          IntList node_of_elem;
          node_of_elem.vide();
          for (int nod=0; nod<nb_som_elem; nod++)
            {
              int nod_glob = elems(num_elem, nod);
              node_of_elem.add_if_not(nod_glob);
            }
          for (int fac=0; fac<nb_faces_elem; fac++)
            {
              int num_fac = elem_face(num_elem, fac);
              // Cerr << "face "<<fac<<" : num_face = "<<num_fac<<finl;
              for (int voisin=0; voisin<2; voisin++) // Two neighbors only
                {
                  int num_elem_v = face_voisins(num_fac,voisin);
                  if ( (num_elem_v!=-1) && (num_elem_v!=num_elem) )
                    {
                      if (aireArray(num_elem_v) > 0. || all_elem_vois) elem_voisins[num_elem].add_if_not(num_elem_v); // neighbor of num_elem via a face
                      // Cerr << "voisin "<<voisin<<" : num_elem_v = "<<num_elem_v<<finl;
                      // Cerr<<"list elem_voisins => ";
                      // for (int v=0; v<elem_voisins[num_elem].size(); v++) Cerr<<(elem_voisins[num_elem])[v]<< " ";
                      // Cerr<<finl;
 
                      // element num_elem_v_v, neighbor of num_elem_v through a face
                      for (int fac_v=0; fac_v<nb_faces_elem; fac_v++)
                        {
                          int num_fac_v = elem_face(num_elem_v, fac_v);
                          for (int voisin_v=0; voisin_v<2; voisin_v++)
                            {
                              int num_elem_v_v = face_voisins(num_fac_v,voisin_v);
                              if ( (num_elem_v_v!=-1) && (num_elem_v_v!=num_elem_v) && (num_elem_v_v!=num_elem) )
                                {
                                  if (aireArray(num_elem_v_v) > 0. || all_elem_vois)
                                    {
                                      int iok = 0 ;
                                      for (int nod_v_v=0; nod_v_v<nb_som_elem; nod_v_v++)
                                        {
                                          int num_nod_v_v = elems(num_elem_v_v, nod_v_v);
                                          if (node_of_elem.contient(num_nod_v_v)) iok = 1;
                                        }
                                      // Cerr << "num_elem_v_v = "<<num_elem_v_v<<" iok = "<<iok<<finl;
                                      if (iok == 1) elem_voisins[num_elem].add_if_not(num_elem_v_v); // neighbor of num_elem via a node
                                    }
                                }
                            }
                        }
                    }
                }
            }
          // Cerr<<"=> "<<elem_voisins[num_elem].size()<<" voisins de "<<num_elem<<" : ";
          // for (int voisin=0; voisin<elem_voisins[num_elem].size(); voisin++) Cerr<<(elem_voisins[num_elem])[voisin]<< " ";
          // Cerr<<finl;
        }
    }
  return;
}
 
double Source_PDF_base::aire_geometrique_IBM(DoubleTab& rotation, int elem)
{
  int dim_esp = Objet_U::dimension;
  const Probleme_base& pb = equation().probleme();
  const Domaine_dis_base& le_dom_dis = pb.domaine_dis();
  const Domaine_VF& the_dom_VF = ref_cast(Domaine_VF,le_dom_dis);
  const Domaine& dom = le_dom_dis.domaine();
  const IntTab& elems = dom.les_elems() ;
  const DoubleTab coordsDom3D=dom.les_sommets();
  int nb_som_elem = dom.nb_som_elem();
 
  DoubleTrav hmax(dim_esp);
  for (int k=0; k < dim_esp; k++)
    {
      double xmin = 1.e30;
      double xmax = -1.e30;
      for (int id_loc=0; id_loc < nb_som_elem; id_loc++)
        {
          int id_glob = elems(elem,id_loc);
          xmin = std::min(xmin, coordsDom3D(id_glob, k));
          xmax = std::max(xmax, coordsDom3D(id_glob, k));
        }
      hmax(k) = xmax-xmin;
    }
  double hcharac = 0.;
  for (int k=0; k<dim_esp; k++) hcharac += hmax(k) * abs(rotation(elem,3*k+(dim_esp-1)));
  hcharac = abs(hcharac);
  // hcharac = (prepro_lu_->get_h_max_elem())(e);
  double aire = the_dom_VF.volumes(elem) / hcharac;
 
  // Cerr<< "=> element: "<<elem<<" ; hmax "<< hmax<<finl;
  // Cerr<< "normal "<<rotation(elem,3*0+(dim_esp-1))<<" "<<rotation(elem,3*1+(dim_esp-1))<<" "<<rotation(elem,3*2+(dim_esp-1))<<finl;
  // Cerr<< "hcharac "<< hcharac<<" ; aire = "<<aire<<finl;
 
  return aire;
}
 
void Source_PDF_base::update_pseudo_level_set_IBM(DoubleTab& vecteur_deplacement, double alpha)
{
  // Update from the pseudo level set and vertex displacement vector
 
  if (!(getInterpolationBool() == true))
    {
      Cerr << "Source_PDF_base::update_pseudo_level_set_IBM: No Interpolation. Aborting..." << finl;
      abort();
    }
  if (!prepro_lu_)
    {
      Cerr << "Source_PDF_base::update_pseudo_level_set_IBM: No IBM prepro. Aborting..." << finl;
      abort();
    }
 
  int dim_esp = Objet_U::dimension;
  assert(dim_esp == 3);
  Interpolation_IBM_base& my_interp = ref_cast_non_const(Interpolation_IBM_base, getInterpolationLu());
  DoubleTab& level_set = ref_cast_non_const(DoubleTab, my_interp.get_pseudo_level_set());
  int nb_som = level_set.dimension(0);
 
  DoubleTab& bary = champ_barycentre_->valeurs();
  DoubleTab& rotation = champ_rotation_->valeurs();
  DoubleTab& aire = champ_aire_->valeurs();
  int nb_elem = aire.dimension(0);
 
  const DoubleTab& nor = my_interp.get_normal_proj_solid();
  assert (nor.dimension(0) == nb_som);
  assert (nor.dimension(1) == dim_esp);
 
  assert (vecteur_deplacement.dimension(0) == nb_som);
  assert (vecteur_deplacement.dimension(1) == dim_esp);
 
  const Probleme_base& pb = equation().probleme();
  const Domaine_dis_base& le_dom_dis = pb.domaine_dis();
  const Domaine& dom = le_dom_dis.domaine();
  const IntTab& elems= dom.les_elems() ;
  int nb_som_elem=dom.nb_som_elem();
  const DoubleTab coordsDom3D=dom.les_sommets();
 
  // Initialization of the level-set
  double nor_nor = 0.;
  for (int e=0; e<nb_elem; e++)
    {
      int nb_nor_nul = 0;
      for (int il=0; il<nb_som_elem; il++)
        {
          int i = elems(e,il);
          if (level_set(i) == -99999.) // level_set not initialized
            {
              nb_nor_nul += 1;
              continue;
            }
          nor_nor = 0.;
          for (int k=0; k<dim_esp; k++) nor_nor += nor(i, k)*nor(i, k);
          if (nor_nor <= 1.e-10) nb_nor_nul += 1; // vertex not computed by PDF
        }
      if ( nb_nor_nul == nb_som_elem )
        for (int il=0; il<nb_som_elem; il++) level_set(elems(e,il)) = -99999.;
    }
 
  // Gradient of phi per element (explicit scheme) + vertex (robustness)
  DoubleTrav grad_phi_e_0(nb_elem,dim_esp) ;
  DoubleTrav grad_phi_i_0(nb_som,dim_esp) ;
  IntTrav contrib_i(nb_som);
  int etarg = -1; // for debug
  // grad_phi per element (normalized)
  for (int e=0; e<nb_elem; e++)
    {
      int initialised_vertex_in_e = 1;
      if (e == etarg) Cerr<<">>> elem = "<<e<<" : grad phi"<<finl;
      // Add contributions within the element
      for (int il=0; il<nb_som_elem; il++)
        {
          int i = elems(e,il);
          if (level_set(i) == -99999.) // level_set not initialized => no contribution
            {
              initialised_vertex_in_e = 0;
              continue;
            }
          int dir_grad_phi_0 = (level_set(i)<0.?-1:(level_set(i)==0.?0.:1)); // grad_phi along increasing phi direction (dir_grad_phi * nor)
          nor_nor = 0.;
          for (int k=0; k<dim_esp; k++) nor_nor += nor(i, k)*nor(i, k);
          if (e == etarg) Cerr<<"i = "<<i<<" norme carree  normal = "<<nor_nor<<" dir_grad_phi_0 = "<<dir_grad_phi_0<<finl;
          if (nor_nor > 1.e-10)
            {
              for (int k=0; k<dim_esp; k++) grad_phi_e_0(e,k) += dir_grad_phi_0*nor(i, k);
            }
        }
      // normalization of grad_phi_e_0 per element if non-null
      // note: if level_set not initialized for all nodes of the element => grad_phi_e_0 is null
      // Treatment of initialised_vertex_in_e = 0 and vertex assembly
      nor_nor = 0.;
      for (int k=0; k<dim_esp; k++) nor_nor += grad_phi_e_0(e,k)*grad_phi_e_0(e,k);
      if (nor_nor > 1.e-10)
        {
          for (int k=0; k<dim_esp; k++) grad_phi_e_0(e,k) /= sqrt(nor_nor);
 
          // If initialised_vertex_in_e = 0
          //   => initialize these vertices from grad_phi_e_0 and the other vertices
          if ( initialised_vertex_in_e == 0)
            {
              // loop to find reference vertex
              double level_set_ref = -99999.;
              DoubleTrav XYZ_ref(dim_esp);
              for (int il=0; il<nb_som_elem; il++)
                {
                  int i = elems(e,il);
                  nor_nor = 0.;
                  for (int k=0; k<dim_esp; k++) nor_nor += nor(i, k)*nor(i, k);
                  if (level_set(i) != -99999. && nor_nor > 1.e-10)
                    {
                      level_set_ref = level_set(i) ;
                      for (int k=0; k<dim_esp; k++) XYZ_ref(k) = coordsDom3D(i,k);
                      break; // found
                    }
                }
              if (e == etarg) Cerr<<"Vertex not initialized : level_set_ref = "<<level_set_ref<<finl;
              // loop to define level_set for uninitialized vertices
              for (int il=0; il<nb_som_elem; il++)
                {
                  int i = elems(e,il);
                  if (level_set(i) == -99999.) // initialize
                    {
                      double prod_sca = 0.;
                      for (int k=0; k<dim_esp; k++) prod_sca += (coordsDom3D(i,k)-XYZ_ref(k))*grad_phi_e_0(e,k);
                      level_set(i) = level_set_ref + prod_sca;
                    }
                }
            }
 
          // Vertex assembly
          for (int il=0; il<nb_som_elem; il++)
            {
              int i = elems(e,il);
              for (int k=0; k<dim_esp; k++) grad_phi_i_0(i,k) += grad_phi_e_0(e,k);
              contrib_i(i) += 1;
            }
        }
 
      if (e == etarg)
        {
          for (int il=0; il<nb_som_elem; il++)
            {
              int i = elems(e,il);
              Cerr<<"i = "<<i<<" ; old level set = "<<level_set(i)<<" ; nor = ";
              for (int k=0; k<dim_esp; k++) Cerr<<nor(i,k)<<" ";
              Cerr<<" ; deplacement = ";
              for (int k=0; k<dim_esp; k++) Cerr<<vecteur_deplacement(i,k)<<" ";
              Cerr<<finl;
            }
          Cerr<<"grad_phi_e_0 = ";
          for (int k=0; k<dim_esp; k++) Cerr<<grad_phi_e_0(e,k)<<" ";
          Cerr<<finl;
        }
    }
 
  // grad_phi per vertex (normalized)
  for (int i=0; i<nb_som; i++)
    {
      if (contrib_i(i) != 0)
        {
          for (int k=0; k<dim_esp; k++) grad_phi_i_0(i,k) /=contrib_i(i);
 
          // normalisation
          nor_nor = 0.;
          for (int k=0; k<dim_esp; k++) nor_nor += grad_phi_i_0(i,k)*grad_phi_i_0(i,k);
          if (nor_nor > 1.e-10)
            for (int k=0; k<dim_esp; k++) grad_phi_i_0(i,k) /= sqrt(nor_nor);
        }
    }
 
  // Advect phi per vertex for all vertices + new_bary for elements crossed by level set zero (arithmetic mean)
  // normal equals grad_phi_e_0 for these elements
  IntTrav advect_faite(nb_som) ;
  aire = 0.; // Zero out the IB area field
  bary = 0.; // Zero out the IB barycenter field
  DoubleTrav new_normal_e(bary);
  for (int e=0; e<nb_elem; e++)
    {
      if (e == etarg) Cerr<<">>> elem = "<<e<<" : Advecte phi"<<finl;
      double nor_nor_e_0 = 0.;
      for (int k=0; k<dim_esp; k++) nor_nor_e_0 += grad_phi_e_0(e,k)*grad_phi_e_0(e,k);
      if (nor_nor_e_0 <= 1.0e-10) continue; // element with undefined element grad phi
 
      double contrib_e = 0.;
      for (int il=0; il<nb_som_elem; il++)
        {
          int i = elems(e,il);
          // Advect phi for vertex i of the element (if not already done)
          if (advect_faite(i) == 0)
            {
              if (e == etarg)
                {
                  Cerr<<" ** node = "<<i<<finl;
                  Cerr<<" vecteur_deplacement = ";
                  for (int k=0; k<dim_esp; k++) Cerr<<vecteur_deplacement(i,k)<<" ";
                  Cerr<<finl;
                }
 
              if (level_set(i) != -99999.)
                {
                  double vdotnorm = 0.0 ;
                  for (int k=0; k<dim_esp; k++) vdotnorm += vecteur_deplacement(i,k)*grad_phi_i_0(i,k);
                  level_set(i) -= alpha*vdotnorm;
                }
              advect_faite(i) = (level_set(i)<0.?-1:1);
 
 
              if (e == etarg)
                {
                  Cerr<<" nor = ";
                  for (int k=0; k<dim_esp; k++) Cerr<<nor(i,k)<<" ";
                  Cerr<<finl;
                  Cerr<<" new level set : "<<level_set(i)<<" ; advect_faite = "<<advect_faite(i)<<finl;
                  Cerr<<" grad_phi_i_0 = ";
                  for (int k=0; k<dim_esp; k++) Cerr<<grad_phi_i_0(i,k)<<" ";
                  Cerr<<finl;
                }
            }
 
          // compute nodal solid projection and barycenter (arithmetic mean)
          for (int k=0; k<dim_esp; k++)
            bary(e,k) += coordsDom3D(i,k)-(level_set(i)*grad_phi_i_0(i,k)); // OP^1 = OX - phi^1 grad_phi^0
          contrib_e += 1.;
          if (e == etarg)
            {
              Cerr<<" ** node = "<<i;
              Cerr<<" ; new level set = "<<level_set(i);
              Cerr<<" ; coords = ";
              for (int k=0; k<dim_esp; k++) Cerr<<coordsDom3D(i,k)<<" ";
              Cerr<<finl;
              Cerr<<" Projection : ";
              for (int k=0; k<dim_esp; k++) Cerr<<coordsDom3D(i,k)-(level_set(i)*grad_phi_i_0(i,k))<<" ";
              Cerr<<finl;
            }
        }
      if (contrib_e != 0.)
        for (int k=0; k<dim_esp; k++) bary(e,k) /= contrib_e;
 
      // If element crossed by level set zero + element grad_phi_e OK => contribution to area field
      int iok = 0;
      // element with level set of different sign?
      double level_ref = 0.;
      for (int il=0; il<nb_som_elem; il++)
        {
          int i = elems(e,il);
          if (level_set(i) != 0. && level_set(i) != -99999.)
            {
              level_ref = level_set(i); // reference: first vertex with phi<>0 and initialized
              break;
            }
        }
      if (level_ref != 0.)
        {
          for (int il=0; il<nb_som_elem; il++)
            {
              int i = elems(e,il);
              if (level_set(i) != 0. && level_set(i) != -99999.)
                {
                  iok = (level_ref*level_set(i) < 0. ? 1:0) ;
                  if (iok == 1)  break;
                }
            }
        }
      // grad_phi_e element OK ?
      if (nor_nor_e_0 > 1.0e-10) iok += 1;
      if (iok == 2)
        {
          // Area definition to be done once the rotations have been updated
          aire(e) = 9999.;
          // Definition of the element normal
          for (int k=0; k<dim_esp; k++) new_normal_e(e, k) = grad_phi_e_0(e,k);
          if (e == etarg)
            {
              Cerr<<"aire = "<<aire(e)<<finl;
              Cerr<<"new_bary_e = ";
              for (int k=0; k<dim_esp; k++) Cerr<<bary(e,k)<<" ";
              Cerr<<finl;
              Cerr<<"new_normal_e = ";
              for (int k=0; k<dim_esp; k++) Cerr<<new_normal_e(e,k)<<" ";
              Cerr<<finl;
            }
        }
      else
        for (int k=0; k<dim_esp; k++) bary(e,k) =0.;
 
    }
  level_set.echange_espace_virtuel();
 
  // Update rotation field per element
  rotation = 0.; // Zero out the rotations
  rotation.echange_espace_virtuel();
  DoubleTrav t1EulerArr(bary);
  DoubleTrav t2EulerArr(bary);
  for (int e=0; e<nb_elem; e++)
    {
      if(aire(e) == 9999.)
        {
          prepro_lu_->computeLocalFrame(new_normal_e, t1EulerArr, t2EulerArr, e);
          prepro_lu_->computeMatRot(new_normal_e, t1EulerArr, t2EulerArr, e);
        }
    }
  rotation = prepro_lu_->get_champ_rotation();
 
  // Rotations up to date => possible modification of the area field
  for (int e=0; e<nb_elem; e++)
    {
      if(aire(e) == 9999.)
        {
          aire(e) = aire_geometrique_IBM(rotation, e);
          if (e == etarg) Cerr<<"aire = "<<aire(e)<<finl;
        }
    }
 
  aire.echange_espace_virtuel();
  bary.echange_espace_virtuel();
  rotation.echange_espace_virtuel();
 
  // finalization: update of the IB
  set_fields_to_prepro(prepro_lu_);
  prepro_lu_->compute_solid_fluid(1);
 
  my_interp.set_fields_from_prepro_to_interp(prepro_lu_);
  my_interp.calculer_normal_et_distance_proj_solid();
  my_interp.set_pseudo_level_set(level_set);
  // my_interp.definir_pseudo_level_set(); //re-initialisation phi
 
  set_variable_imposee();
  prepro_lu_-> Save_Med_File();
 
}
 
void Source_PDF_base::update_elem_IBM(DoubleTab& vecteur_deplacement, double alpha, double raid)
{
  // Update from the elemental barycenters of Source_PDF_base
 
  int dim_esp = Objet_U::dimension;
  assert (dim_esp == 3);
 
  DoubleTab& aire = ref_cast_non_const(DoubleTab, champ_aire_->valeurs());
  DoubleTab& bary = ref_cast_non_const(DoubleTab, champ_barycentre_->valeurs());
  DoubleTab& rotation = ref_cast_non_const(DoubleTab, champ_rotation_->valeurs());
  const Probleme_base& pb = equation().probleme();
  const Domaine_dis_base& le_dom_dis = pb.domaine_dis();
  const Domaine& dom = le_dom_dis.domaine();
  int nb_elem_tot = dom.nb_elem_tot();
  const DoubleTab coordsDom3D=dom.les_sommets();
 
  int nb_elem = aire.dimension(0);
  assert (nb_elem == vecteur_deplacement.dimension(0));
  assert (dim_esp == vecteur_deplacement.dimension(1));
  IntTab indic_dead_cell(nb_elem);
  indic_dead_cell = 0;
 
  // compute neighbors of each crossed element
  IntLists elem_voisins(nb_elem_tot);
  compute_NeighNode_IBM_elem(aire, elem_voisins);
 
  // update the barycenter field
  IntTrav toward_new_elem(nb_elem);
 
  // stpe1 : prediction
  DoubleTrav bary_zero(bary);
  bary_zero = bary;
  double norm_vd2;
  for (int e=0; e<nb_elem; e++)
    {
      norm_vd2 = 0.;
      for (int k=0; k<dim_esp; k++) norm_vd2 += vecteur_deplacement(e,k)*vecteur_deplacement(e,k);
      if (aire(e)>0.)
        {
          for (int k=0; k<dim_esp; k++) bary(e,k) = bary_zero(e,k) + alpha * vecteur_deplacement(e,k) ;
          if (elem_voisins[e].size() >= 2 && sqrt(norm_vd2) > 1.e-10)
            {
              for (int voi=0; voi<elem_voisins[e].size(); voi++)
                {
                  for (int k=0; k<dim_esp; k++) bary(e,k) += alpha * raid * bary_zero((elem_voisins[e])[voi],k);
                }
              for (int k=0; k<dim_esp; k++) bary(e,k) /= (1 + alpha * raid * elem_voisins[e].size());
            }
          toward_new_elem(e) = dom.chercher_elements(bary(e,0),bary(e,1),bary(e,2));
          if (aire(toward_new_elem(e)) <= 0. ) indic_dead_cell(toward_new_elem(e)) = 1;
        }
      else //to avoid null being considered as element number 0
        {
          toward_new_elem(e) =-2; // rien a deplacer
        }
    }
  indicateur_dead_cell_ = indic_dead_cell;
 
  // step2 : displacement of fields towards new elements
  DoubleTrav new_aire(aire);
  DoubleTrav new_bary(bary);
  DoubleTrav new_rotation(rotation);
  IntTrav contrib(nb_elem);
 
  for (int e=0; e<nb_elem; e++)
    {
      if (toward_new_elem(e) == -1 ) //element non trouve
        {
          Cerr<<"Source_PDF_base::update_elem_IBM: element displacement is not valid."<<finl;
          exit();
        }
      else if (toward_new_elem(e) == -2 ) // element not crossed by gamma
        {
        }
      else //displacement of data (barycenters and rotation matrices)
        {
          for (int k=0; k<dim_esp; k++) new_bary(toward_new_elem(e),k) += bary(e,k);
          for (int k=0; k<rotation.dimension(1); k++) new_rotation(toward_new_elem(e),k) += rotation(e,k);
          new_aire(toward_new_elem(e)) += aire(e);
          contrib(toward_new_elem(e)) += 1;
        }
    }
 
  // Arithmetic mean for barycenters and rotation matrices; arithmetic mean or geometric definition for area
 
  int aire_geometriq_ok = 1;
 
  for (int e=0; e<nb_elem; e++)
    {
      if(contrib(e)>0)
        {
          for (int k=0; k<dim_esp; k++) new_bary(e,k) /= contrib(e);
          for (int k=0; k<rotation.dimension(1); k++) new_rotation(e,k) /= contrib(e);
 
          if (aire_geometriq_ok)
            new_aire(e) = aire_geometrique_IBM(new_rotation, e);
          else new_aire(e) /= contrib(e);
 
        }
    }
 
  // Maj champs
  for (int e=0; e<nb_elem; e++)
    {
      aire(e) = new_aire(e);
      for (int k=0; k<dim_esp; k++) bary(e,k) = new_bary(e,k);
      for (int k=0; k<rotation.dimension(1); k++) rotation(e,k) = new_rotation(e,k);
    }
 
  // Verification of topology conservation
  // Treatment of elements having 0 or 1 neighbor
  // the neighbors of the predecessor of elem must be neighbors of elem
 
  // compute neighbors of each new crossed element
  IntLists elem_voisins_new(nb_elem_tot);
  compute_NeighNode_IBM_elem(aire, elem_voisins_new);
 
  // find displaced elements with fewer than 2 neighbors
  for (int e=0; e<nb_elem_tot; e++)
    {
      if ( elem_voisins_new[e].size() < 2)
        {
          IntList antecedents;
          antecedents.vide();
          // list of elements (antecedents) that have moved towards e
          for (int e_old=0; e_old<nb_elem_tot; e_old++)
            if (toward_new_elem(e_old) == e) antecedents.add_if_not(e_old);
//          if (antecedents.size() != 0) Cerr<<"element "<<e<<" => "<<finl;
          IntList to_link;
          to_link.vide();
          // loop over predecessors of e
          for (int ant=0; ant<antecedents.size(); ant++)
            {
              // Cerr<<" voisins de antecedent = "<<antecedents[ant] <<" ";
              // loop over old neighbors of each predecessor of e
              int nb_vois_ant = elem_voisins[antecedents[ant]].size();
              for (int vois_ant=0; vois_ant<nb_vois_ant; vois_ant++)
                {
                  int old_vois = (elem_voisins[antecedents[ant]])[vois_ant];
                  // is the old neighbor still a neighbor of e?
                  int i_present = 0;
                  for (int vois_e=0; vois_e<elem_voisins_new[e].size(); vois_e++)
                    if ( (elem_voisins_new[e])[vois_e] == toward_new_elem(old_vois) ) i_present = 1;
                  if ( (!i_present) && (toward_new_elem(old_vois) != e) ) to_link.add_if_not(toward_new_elem(old_vois));
                  // Cerr<<old_vois<<" ";
                }
              // Cerr<<finl;
 
              // if (to_link.size() != 0)  Cerr<<"to_link ";
              DoubleTrav Xtest(dim_esp);
              for (int n=0; n <to_link.size(); n++)
                {
                  // Cerr<<to_link[n]<<" ";
                  // Xdeb (bary(e,0),bary(e,1),bary(e,2));
                  // Xfin (bary(to_link[n],0),bary(to_link[n],1),bary(to_link[n],2));
                  int decoup = 5;
                  for (int nb_step=0; nb_step<(decoup-1); nb_step++)
                    {
                      for (int k=0; k<dim_esp; k++) Xtest(k) = bary(e,k) + (bary(to_link[n], k) - bary(e,k))/decoup*(nb_step+1);
                      int elem_Xtest = dom.chercher_elements(Xtest(0),Xtest(1),Xtest(2));
                      if (elem_Xtest != -1)
                        {
                          if (aire(elem_Xtest) <= 0.)
                            {
                              for (int k=0; k<dim_esp; k++) bary(elem_Xtest,k) = Xtest(k);
                              for (int k=0; k<rotation.dimension(1); k++) rotation(elem_Xtest,k) = (rotation(e,k) + rotation(to_link[n],k)) / 2.;
                              aire(elem_Xtest) = aire_geometrique_IBM(rotation, elem_Xtest);
                            }
                        }
                    }
                }
              // Cerr<<finl;
            }
        }
 
    }
 
  // finalisation
  aire.echange_espace_virtuel();
  bary.echange_espace_virtuel();
  rotation.echange_espace_virtuel();
 
 
  if (prepro_lu_)
    {
      set_fields_to_prepro(prepro_lu_);
      prepro_lu_->compute_solid_fluid(1);
      if(getInterpolationBool() == true)
        {
          Interpolation_IBM_base& my_interp = ref_cast_non_const(Interpolation_IBM_base, getInterpolationLu());
          my_interp.set_fields_from_prepro_to_interp(prepro_lu_);
        }
    }
  set_variable_imposee();
}
 
void Source_PDF_base::verify_results_prepro()
{
  if (!prepro_lu_) return;
  Cerr<<"(IBM) source PDF base: Prepro IBM verification if any ...."<<finl;
 
  const DoubleTab& aireArray = prepro_lu_->champ_aire_->valeurs();
  const DoubleTab& baryArray = prepro_lu_->champ_bary_->valeurs();
  const DoubleTab& normalArray = prepro_lu_->champ_normal_->valeurs();
  const DoubleTab& rotationArray = prepro_lu_->champ_rotation_->valeurs();
  const DoubleTab& isNodeDirichletArray = prepro_lu_->isNodeDirichlet_->valeurs();
  const DoubleTab& solideArray = prepro_lu_->solid_points_->valeurs();
  const DoubleTab& fluidArray = prepro_lu_->fluid_points_->valeurs();
 
  const DoubleTab& solid_elemsArray = prepro_lu_->solid_elems_->valeurs();
  const DoubleTab& fluid_elemsArray = prepro_lu_->fluid_elems_->valeurs();
  const DoubleTab& corresp_elemsArray = prepro_lu_->corresp_elems_->valeurs();
 
  int dim_esp = Objet_U::dimension;
 
  ///////////////////////////////// Base fields /////////////////////////////////
  // Area
  if (champ_aire_lu_)
    {
      const DoubleTab& aireArray_src_pdf = champ_aire_lu_->valeurs();
      assert(aireArray_src_pdf.dimension(0)==aireArray.dimension(0));
      Cerr<<"Prepro_IBM::verify_results_prepro: relative L2-norm and abs max-norm errors for aire :"<<finl;
      comp_diff_L2_max(aireArray.dimension(0), aireArray.dimension(1), aireArray_src_pdf, aireArray);
    }
 
  // Barycentre
  if (champ_barycentre_lu_)
    {
      const DoubleTab& baryArray_src_pdf = champ_barycentre_lu_->valeurs();
      assert(baryArray_src_pdf.dimension(0)==baryArray.dimension(0));
      assert(baryArray_src_pdf.dimension(1)==dim_esp);
      Cerr<<"Prepro_IBM::verify_results_prepro: relative L2-norm and abs max-norm errors for barycenter :"<<finl;
      comp_diff_L2_max(baryArray.dimension(0), dim_esp, baryArray_src_pdf, baryArray);
    }
 
  // Normale
  if (champ_rotation_lu_)
    {
      const DoubleTab& rotArray_src_pdf = champ_rotation_lu_->valeurs();
      assert(rotArray_src_pdf.dimension(1)==dim_esp*dim_esp);
      assert(rotArray_src_pdf.dimension(0)==normalArray.dimension(0));
      assert(normalArray.dimension(1)==dim_esp);
      Cerr<<"Prepro_IBM::verify_results_prepro: relative L2-norm and abs max-norm errors for normal :"<<finl;
      for (int k = 0; k < dim_esp; k++)
        {
          double err_L2_norm =0.;
          double err_max_norm = 0.;
          double L2_norm = 0.;
          for (int elem = 0; elem <normalArray.dimension(0); elem++)
            if (aireArray(elem) > 0.)
              {
                double diff_norm = abs(normalArray(elem,k) - rotArray_src_pdf(elem,3*k+(dim_esp-1)));
                err_L2_norm  += diff_norm * diff_norm ;
                L2_norm += rotArray_src_pdf(elem,3*k+(dim_esp-1)) * rotArray_src_pdf(elem,3*k+(dim_esp-1)) ;
                if (diff_norm > err_max_norm) err_max_norm = diff_norm;
              }
          if (L2_norm > 1.0e-12) err_L2_norm = sqrt(err_L2_norm / L2_norm);
 
          Cerr<<"    composant # "<<k<<" => "<<err_L2_norm<<" "<<err_max_norm;
        }
      Cerr<<finl;
    }
 
  // Rotation
  if (champ_rotation_lu_)
    {
      const DoubleTab& rotArray_src_pdf = champ_rotation_lu_->valeurs();
      assert(rotArray_src_pdf.dimension(1)==dim_esp*dim_esp);
      assert(rotArray_src_pdf.dimension(0)==rotationArray.dimension(0));
      assert(rotationArray.dimension(1)==dim_esp*dim_esp);
      Cerr<<"Prepro_IBM::verify_results_prepro: relative L2-norm and abs max-norm errors for rotation :"<<finl;
      comp_diff_L2_max(rotationArray.dimension(0), dim_esp*dim_esp, rotArray_src_pdf, rotationArray);
    }
 
  ///////////////////////////////// Interpolation /////////////////////////////////
 
  if(interpolation_bool_)
    {
      Cerr<<"Prepro_IBM::verify_results_prepro: Interpolation_IBM"<<finl;
 
      // is_dirichlet
      if (interpolation_lue_->is_dirichlet_lu_)
        {
          const DoubleTab& isNodeDirichletArray_src_pdf = interpolation_lue_->is_dirichlet_lu_->valeurs();
          assert(isNodeDirichletArray_src_pdf.dimension(0)==isNodeDirichletArray.dimension(0));
          assert(isNodeDirichletArray_src_pdf.dimension(1)==isNodeDirichletArray.dimension(1));
          Cerr<<"Prepro_IBM::verify_results_prepro: relative L2-norm and abs max-norm errors for isNodeDirichlet :"<<finl;
          comp_diff_L2_max(isNodeDirichletArray.dimension(0), isNodeDirichletArray.dimension(1), isNodeDirichletArray_src_pdf, isNodeDirichletArray);
        }
 
      // Projection solide
      if (interpolation_lue_->solid_points_lu_)
        {
          const DoubleTab& solideArray_src_pdf = interpolation_lue_->solid_points_lu_->valeurs();
          assert(solideArray_src_pdf.dimension(0)==solideArray.dimension(0));
          assert(solideArray_src_pdf.dimension(1)==dim_esp);
          Cerr<<"Prepro_IBM::verify_results_prepro: relative L2-norm and abs max-norm errors for solid projection :"<<finl;
          comp_diff_L2_max(solideArray.dimension(0), dim_esp, solideArray_src_pdf, solideArray);
        }
 
      // correspondance elems
      if (interpolation_lue_->corresp_elems_lu_)
        {
          const DoubleTab& corresp_elems_src_pdf = interpolation_lue_->corresp_elems_lu_->valeurs();
          assert(corresp_elems_src_pdf.dimension(1)==1);
          assert(corresp_elems_src_pdf.dimension(0)==corresp_elemsArray.dimension(0));
          assert(corresp_elemsArray.dimension(1)==1);
          Cerr<<"Prepro_IBM::verify_results_prepro: relative L2-norm and abs max-norm errors for correspondance elems :"<<finl;
          comp_diff_L2_max(corresp_elemsArray.dimension(0), 1, corresp_elems_src_pdf, corresp_elemsArray);
        }
 
      Cerr<<"interpolation_lue_->que_suis_je() = "<<interpolation_lue_->que_suis_je()<<finl;
      // Interpolation_IBM_mean_gradient
      if (interpolation_lue_->que_suis_je() == "Interpolation_IBM_gradient_moyen")
        {
          Cerr<<"Prepro_IBM::verify_results_prepro: >> mean_gradient"<<finl;
          Interpolation_IBM_mean_gradient& interp = ref_cast(Interpolation_IBM_mean_gradient,interpolation_lue_.valeur());
 
          // Element projection solide
          if (interp.solid_elems_lu_)
            {
              const DoubleTab& solid_elemsArray_src_pdf = interp.solid_elems_lu_->valeurs();
              assert(solid_elemsArray_src_pdf.dimension(0)==solid_elemsArray.dimension(0));
              assert(solid_elemsArray_src_pdf.dimension(1)==solid_elemsArray.dimension(1));
              Cerr<<"Prepro_IBM::verify_results_prepro: relative L2-norm and abs max-norm errors for solid projection element :"<<finl;
              comp_diff_L2_max(solid_elemsArray.dimension(0), solid_elemsArray.dimension(1), solid_elemsArray_src_pdf, solid_elemsArray);
            }
        }
 
      // Interpolation_IBM_elem_fluid
      if (interpolation_lue_->que_suis_je() == "Interpolation_IBM_element_fluide")
        {
          Cerr<<"Prepro_IBM::verify_results_prepro: >> elem_fluide"<<finl;
          Interpolation_IBM_elem_fluid& interp = ref_cast(Interpolation_IBM_elem_fluid,interpolation_lue_.valeur());
 
          // Projection fluide
          if (interp.fluid_points_lu_)
            {
              const DoubleTab& fluidArray_src_pdf = interp.fluid_points_lu_->valeurs();
              assert(fluidArray_src_pdf.dimension(0)==fluidArray.dimension(0));
              assert(fluidArray_src_pdf.dimension(1)==dim_esp);
              Cerr<<"Prepro_IBM::verify_results_prepro: relative L2-norm and abs max-norm errors for fluid projection :"<<finl;
              comp_diff_L2_max(fluidArray.dimension(0), dim_esp, fluidArray_src_pdf, fluidArray);
            }
 
          // Element projection fluide
          if (interp.fluid_elems_lu_)
            {
              const DoubleTab& fluid_elemsArray_src_pdf = interp.fluid_elems_lu_->valeurs();
              assert(fluid_elemsArray_src_pdf.dimension(0)==fluid_elemsArray.dimension(0));
              assert(fluid_elemsArray_src_pdf.dimension(1)==fluid_elemsArray.dimension(1));
              Cerr<<"Prepro_IBM::verify_results_prepro: relative L2-norm and abs max-norm errors for fluid projection element :"<<finl;
              comp_diff_L2_max(fluid_elemsArray.dimension(0), fluid_elemsArray.dimension(1), fluid_elemsArray_src_pdf, fluid_elemsArray);
            }
        }
 
    }
}
 
void Source_PDF_base::filtre_CLD(DoubleTab& filtre) const
{
  Cerr << "Source_PDF_base: Not implemented for current discretisation. Aborting..." << finl;
  abort();
}
 
void Source_PDF_base::rotate_imposed_velocity(DoubleTab& vitesse_imposee)
{
  const Domaine_dis_base& Domaine_dis = equation().probleme().domaine_dis();
  int nb_elem_tot=Domaine_dis.domaine().nb_elem_tot();
  DoubleTab& rotation = champ_rotation_->valeurs();
 
  // VDF/VEF by default
  const Domaine_VF& the_dom = ref_cast(Domaine_VF,Domaine_dis);
  const IntTab& elems = (equation().discretisation().is_ef() ? Domaine_dis.domaine().les_elems() : the_dom.elem_faces());
  int nb_som_tot = the_dom.nb_faces_tot();
  int nb_dof_elem = Domaine_dis.domaine().nb_faces_elem();
  if (equation().discretisation().is_ef())
    {
      nb_dof_elem = Domaine_dis.domaine().nb_som_elem();
      nb_som_tot= Domaine_dis.domaine().nb_som_tot();
    }
  else if (!(equation().discretisation().is_vdf()) && !(equation().discretisation().is_vef()))
    {
      Cerr<<"Source_PDF_base::rotate_imposed_velocit: discretisation is not EF, VEF, VDF : "<<equation().discretisation()<<finl;
    }
 
  // COPIES
  DoubleTab vitesse = vitesse_imposee;
  assert(Objet_U::dimension==3);
  int dim_var=vitesse_imposee.dimension(1);
  assert(Objet_U::dimension==dim_var);
 
  //vitesse_imposee.echange_espace_virtuel();
  DoubleTrav marqueur=equation().probleme().get_champ("vitesse").valeurs();
  vitesse_imposee = 0.;
  marqueur = 0.;
 
  for (int num_elem=0; num_elem<nb_elem_tot; num_elem++)
    {
      double norm =  sqrt(rotation(num_elem,0)*rotation(num_elem,0)
                          +rotation(num_elem,1)*rotation(num_elem,1)
                          +rotation(num_elem,2)*rotation(num_elem,2)
                          +rotation(num_elem,3)*rotation(num_elem,3)
                          +rotation(num_elem,4)*rotation(num_elem,4)
                          +rotation(num_elem,5)*rotation(num_elem,5)
                          +rotation(num_elem,6)*rotation(num_elem,6)
                          +rotation(num_elem,7)*rotation(num_elem,7)
                          +rotation(num_elem,8)*rotation(num_elem,8));
      if (norm > 1e-6)
        {
          for (int i=0; i<nb_dof_elem; i++)
            {
              int s = elems(num_elem,i);
              for (int c=0; c<dim_var; c++)
                {
                  for (int k=0; k<dim_var; k++)
                    {
                      vitesse_imposee(s,c) += rotation(num_elem,3*c+k)*vitesse(s,k);
                      //Cerr << "k = " << k << ", c = " << c << ", 3*c+k = " << 3*c+k << ", rotation = " << rotation(num_elem,3*c+k) << finl;
                    }
                }
              marqueur(s,0) += 1.0;
            }
        }
    }
  for (int i=0; i<nb_som_tot; i++)
    {
      if(marqueur(i,0) > 0.)
        {
          for (int c=0; c<dim_var; c++)
            {
              vitesse_imposee(i,c) /= marqueur(i,0);
            }
        }
    }
}
 
void Source_PDF_base::compute_variable_imposee_projete(const DoubleTab& marqueur, const DoubleTab& points, double val, double eps)
{
  const Domaine_dis_base& le_dom_dis = equation().probleme().domaine_dis();
 
  int nb_som=le_dom_dis.domaine().nb_som();
  int dim = Objet_U::dimension;
  const DoubleTab& coords = le_dom_dis.domaine().coord_sommets();
  ArrOfDouble x(dim);
  int dim_var = equation().inconnue().valeurs().dimension(1);
  for (int i = 0; i < nb_som; i++)
    {
      if ((marqueur(i) >= val + eps) || (marqueur(i) <= val - eps))
        {
          for (int j = 0; j < dim; j++)
            {
              x[j] = points(i,j);
            }
          for (int j = 0; j < dim_var; j++)
            {
              modele_lu_.variable_imposee_->valeurs()(i,j) = modele_lu_.get_variable_imposee(x,j);
            }
        }
      else
        {
          for (int j = 0; j < dim; j++)
            {
              x[j] = coords(i,j);
            }
          for (int j = 0; j < dim_var; j++)
            {
              modele_lu_.variable_imposee_->valeurs()(i,j) = modele_lu_.get_variable_imposee(x,j);
            }
        }
    }
}
 
void Source_PDF_base::associer_domaines(const Domaine_dis_base& domaine_dis,
                                        const Domaine_Cl_dis_base& domaine_Cl_dis)
{
  Cerr << "Source_PDF_base: Not implemented for current discretisation. Aborting..." << finl;
  abort();
}
 
void Source_PDF_base::compute_indicateur_nodal_champ_aire()
{
  // indicateur_nodal_champ_aire_ = 1 if vertex belongs to an element with non-zero area
  // 0 otherwise
  const DoubleTab& aire=champ_aire_->valeurs();
  const Domaine_dis_base& Domaine_dis = equation().probleme().domaine_dis();
  int nb_elems=Domaine_dis.domaine().nb_elem_tot();
  int nb_nodes=Domaine_dis.domaine().nb_som_tot();
  int nb_som_elem=Domaine_dis.domaine().nb_som_elem();
  const IntTab& elems= Domaine_dis.domaine().les_elems() ;
 
  DoubleTab indic(nb_nodes);
  indic = 0.;
  for (int num_elem=0; num_elem<nb_elems; num_elem++)
    {
      if (aire(num_elem)>0.)
        {
          for (int i=0; i<nb_som_elem; i++)
            {
              int s1=elems(num_elem,i);
              indic(s1) = 1.;
            }
        }
    }
  indicateur_nodal_champ_aire_ = indic;
}
 
double Source_PDF_base::fonct_regul_PDF(const int elem, const double dist)
{
  double coeff = 1.;
 
  int size_dead_cell_tab = indicateur_dead_cell_.size();
  if (size_dead_cell_tab == 0) return coeff;
  if (elem >= size_dead_cell_tab)
    {
      Cerr<<"Source_PDF_base::fonct_regul_PDF: element = "<<elem<<" >= size of tab indicateur_dead_cell = "<<size_dead_cell_tab<<finl;
      exit();
    }
  if (indicateur_dead_cell_(elem) != 1) return coeff;
 
  const Probleme_base& pb = equation().probleme();
  const Domaine_dis_base& le_dom_dis = pb.domaine_dis();
  const Domaine_VF& the_dom_VF = ref_cast(Domaine_VF,le_dom_dis);
  DoubleTab& rotation = ref_cast_non_const(DoubleTab, champ_rotation_->valeurs());
 
  double h_max_elem = the_dom_VF.volumes(elem) / aire_geometrique_IBM(rotation, elem);
  double coeff_regul_PDF = modele_lu_.regul_coeff_PDF_;
  double h_ref = coeff_regul_PDF * h_max_elem;
  coeff = max((h_ref - dist)/h_ref , 0.);
 
  return coeff;
}
 
double Source_PDF_base::fonct_coeff(const double rho_m, const double aire, const double dt) const
{
  // return = 0 if area in element <= 0
  // return = rho/dt * coeff_relax otherwise
  double val_coeff = 0.;
  if (aire<=0.)
    {
      return val_coeff;
    }
 
  double inv_dt = 1./dt;
  // Cerr<<"dt for Source_PDF_base::fonct_coeff : "<< dt <<finl;
  val_coeff = rho_m * inv_dt;
  const double tps_cour = equation().probleme().schema_temps().temps_courant();
  if (temps_relax_ != 1.0e+12)
    {
      // double bco = 100.*(1.-exp(-0.1/echelle_relax_)); // tps_cour_intersection = 0.1 * temps_relax_
      double deltrel = tps_cour - temps_relax_;
      double coeff_relax = (1.+tanh(deltrel/(echelle_relax_ * temps_relax_)))/2.;
      // double coeff_relax = min (1., std::max(1.-bco*deltrel*deltrel/temps_relax_/temps_relax_, exp(deltrel/(echelle_relax_ * temps_relax_))));
      // double coeff_relax = 1.0 - exp(-tps_cour / temps_relax_);
      // double coeff_relax = tanh(tps_cour / temps_relax_);
      val_coeff *= coeff_relax ;
      // Cerr<< "coeff. relax. for PDF = "<<coeff_relax<<" val_coeff = "<<val_coeff <<finl;
    }
  return val_coeff;
}
 
// TENSOR CALCULATION
// To possibly handle space directions differently
ArrOfDouble Source_PDF_base::get_tuvw_local() const
{
  assert(Objet_U::dimension==3);
  ArrOfDouble tuvw(dimension) ;
  tuvw[0] = 1.0 / modele_lu_.eta_;
  tuvw[1] = 1.0 / modele_lu_.eta_;
  tuvw[2] = 1.0 / modele_lu_.eta_;
  return tuvw ;
}
 
DoubleVect Source_PDF_base::diag_coeff_elem(ArrOfDouble& variable_elem, const DoubleTab& rotation, int num_elem) const
{
  assert(Objet_U::dimension==3);
  ArrOfDouble tuvw(dimension);
 
  tuvw = get_tuvw_local();
 
  ArrOfDouble sum_dir_loc(dimension) ;
  for (int i=0; i<dimension; i++)
    {
      sum_dir_loc[i]=0.;
      for (int k=0; k<dimension; k++)
        sum_dir_loc[i]+=rotation(num_elem,3*i+k);
    }
 
  DoubleVect diag_coef(dimension);
  for (int comp=0; comp<dimension; comp++)
    {
      diag_coef(comp) = 0. ;
      for (int k=0; k<dimension; k++)
        diag_coef(comp)+=tuvw[k]*sum_dir_loc[k]*rotation(num_elem,3*k+comp);
    }
 
  return diag_coef ;
}
 
DoubleTab Source_PDF_base::compute_coeff_elem() const
{
  // coeff = 1 if area in element <= 0
  // coeff = 1 + (Ksi/eta) * coeff_relax otherwise
  const Domaine_dis_base& domaine_dis = equation().probleme().domaine_dis();
  const IntTab& elems= domaine_dis.domaine().les_elems() ;
  int nb_som_elem=domaine_dis.domaine().nb_som_elem();
  int nb_elems=domaine_dis.domaine().nb_elem_tot();
  const DoubleTab& rotation=champ_rotation_->valeurs();
  const DoubleTab& aire = champ_aire_->valeurs();
 
  int dim_var = equation().inconnue().valeurs().dimension(1);
  ArrOfDouble variable_elem(dim_var);
  int dim_esp = Objet_U::dimension ;
  if (dim_var != 1 && dim_var != dim_esp)
    {
      Cerr<<"compute_coeff_elem: variable dimension different from 1 or "<<dim_esp<<"! "<<finl;
      exit();
    }
  DoubleTab coeff(nb_elems, dim_var);
 
  double dt = (ref_cast_non_const(Source_PDF_base, *this)).calcul_dt_pdf();
  const DoubleTab& rho_m=champ_rho_->valeurs();
  if (equation().probleme().schema_temps().temps_courant()==0)
    {
      dt = 1.0;
    }
  // Cerr<<"dt for compute_coeff_elem : "<< dt <<finl;
 
  double val_coeff ;
  DoubleVect val_diag_coef(dim_var);
  for (int num_elem=0; num_elem<nb_elems; num_elem++)
    {
      if (aire(num_elem)<=0.)
        {
          for (int comp=0; comp<dim_var; comp++) coeff(num_elem, comp) = 1. ;
        }
      else
        {
          const DoubleTab& variable=equation().inconnue().valeurs();
          variable_elem=0;
          for (int s=0; s<nb_som_elem; s++)
            {
              int som_glob=elems(num_elem,s);
              for (int comp=0; comp<dim_var; comp++)
                variable_elem[comp]+=variable(som_glob,comp);
            }
          // rho/dt * coeff_relax
          val_coeff = fonct_coeff(rho_m(num_elem), aire(num_elem), dt) ;
          // 1/eta
          if (dim_var == dim_esp)
            {
              val_diag_coef = diag_coeff_elem(variable_elem, rotation, num_elem) ;
            }
          else
            {
              val_diag_coef(0) = 1.0 / modele_lu_.eta_ ;
            }
          // coeff = 1 + (Ksi/eta) * coeff_relax
          for (int comp=0; comp<dim_esp; comp++)
            coeff(num_elem, comp) = 1.+ val_coeff * val_diag_coef(comp) * dt / rho_m(num_elem) ;
        }
    }
 
  return coeff;
}
 
DoubleTab Source_PDF_base::compute_coeff_matrice() const
{
  // coeff vertex = Sigma_elem coeff_elem / Nb contributions;  with:
  // coeff_elem = Ksi/eta * coeff_relax if element with area <>= 0
  // coeff_elem = 0 otherwise
  const Domaine_dis_base& Domaine_dis = equation().probleme().domaine_dis();
  int nb_elems=Domaine_dis.domaine().nb_elem_tot();
 
  // VDF/VEF by default
  const Domaine_VF& the_dom = ref_cast(Domaine_VF,Domaine_dis);
  const IntTab& elems = (equation().discretisation().is_ef() ? Domaine_dis.domaine().les_elems() : the_dom.elem_faces());
  int nb_dof_tot = the_dom.nb_faces_tot();
  int nb_dof_elem = Domaine_dis.domaine().nb_faces_elem();
  if (equation().discretisation().is_ef())
    {
      nb_dof_tot= Domaine_dis.domaine().nb_som_tot();
      nb_dof_elem = Domaine_dis.domaine().nb_som_elem();
    }
  else if (!(equation().discretisation().is_vdf()) && !(equation().discretisation().is_vef()))
    {
      Cerr<<"Source_PDF_base::compute_coeff_matrice: discretisation is not EF, VEF, VDF : "<<equation().discretisation()<<finl;
    }
 
  const DoubleTab& rotation=champ_rotation_->valeurs();
  const DoubleTab& aire = champ_aire_->valeurs();
 
  int dim_esp = Objet_U::dimension ;
  const DoubleTab& variable=equation().inconnue().valeurs();
  int dim_var = variable.dimension(1);
  if (dim_var != 1 && dim_var != dim_esp)
    {
      Cerr<<"compute_coeff_matrice: variable dimension different from 1 or "<<dim_esp<<"! "<<finl;
      exit();
    }
  ArrOfDouble variable_elem(dim_var);
  DoubleTab coeff(variable);
  IntTab contrib(nb_dof_tot, dim_var);
  contrib = 0;
  coeff = 0.;
 
  double dt = (ref_cast_non_const(Source_PDF_base, *this)).calcul_dt_pdf();
  const DoubleTab& rho_m=champ_rho_->valeurs();
  // Cerr<<"dt for compute_coeff_matrice : "<< dt <<finl;
 
  double val_coeff ;
  DoubleVect val_diag_coef(dim_var), coeff_el(dim_var) ;
  for (int num_elem=0; num_elem<nb_elems; num_elem++)
    {
      if (aire(num_elem)<=0.)
        {
          for (int comp=0; comp<dim_var; comp++) coeff_el(comp) = 0. ;
        }
      else
        {
          variable_elem=0.;
          for (int s=0; s<nb_dof_elem; s++)
            {
              int dof_glob=elems(num_elem,s);
              for (int comp=0; comp<dim_var; comp++)
                variable_elem[comp]+=variable(dof_glob,comp)/nb_dof_elem;
            }
          // rho/dt * coeff_relax
          val_coeff = fonct_coeff(rho_m(num_elem), aire(num_elem), dt) ;
          //  1/eta
          if (dim_var == dim_esp)
            {
              val_diag_coef = diag_coeff_elem(variable_elem, rotation, num_elem) ;
            }
          else
            {
              val_diag_coef(0) = 1.0 / modele_lu_.eta_ ;
            }
          // coeff = Ksi/eta * coeff_relax
          for (int comp=0; comp<dim_var; comp++)
            coeff_el(comp) = val_coeff * val_diag_coef(comp) * dt / rho_m(num_elem) ;
        }
      for (int i = 0; i<nb_dof_elem; i++)
        {
          int s = elems(num_elem,i);
          for (int comp=0; comp<dim_var; comp++)
            {
              if (coeff_el(comp) > 0.)
                {
                  coeff(s,comp)+= coeff_el(comp);
                  contrib(s, comp) += 1;
                }
              // Cerr<<"Source_PDF_EF: coeff_el ("<<comp<<"): "<<" "<<coeff_el(comp)<<finl;
            }
        }
    }
 
  double val = 0.;
  for (int s = 0; s<nb_dof_tot; s++)
    {
      for (int comp=0; comp<dim_var; comp++)
        {
          val = 1.;
          if (contrib(s,comp) != 0)
            {
              val += coeff(s,comp)/contrib(s,comp);
            }
          coeff(s,comp) = val;
          // Cerr<<"Source_PDF_EF: contrib coeff ("<<s<<","<<comp<<"): "<<contrib(s,comp)<<" "<<coeff(s,comp)<<finl;
        }
    }
  // Cerr << "Coeff matrice" << coeff << finl;
  return coeff;
}
 
void Source_PDF_base::multiply_coeff_volume(DoubleTab& coeff) const
{
  Cerr << "Source_PDF_base: Not implemented for current discretisation. Aborting..." << finl;
  abort();
}
 
DoubleTab Source_PDF_base::compute_pond(const DoubleTab& rho_m, const DoubleTab& aire, const DoubleVect& volume, int& coef1, int& nb_elems) const
{
// return = 0 if area in element <= 0
// return = rho/dt * coeff_relax * volume / (coef1*coef1)  otherwise
  DoubleTab pond = rho_m;
  double inv_coef1 = 1. / coef1;
 
  double dt = (ref_cast_non_const(Source_PDF_base, *this)).calcul_dt_pdf();
 
  for (int num_elem=0; num_elem<nb_elems; num_elem++)
    {
      pond(num_elem) = fonct_coeff(rho_m(num_elem), aire(num_elem), dt) ;
      pond(num_elem) *= volume(num_elem)*inv_coef1*inv_coef1;
    }
 
  return pond;
}
 
DoubleTab& Source_PDF_base::ajouter(DoubleTab& secmem) const
{
  if(has_interface_blocs())
    {
      ajouter_blocs({}, secmem);
      return secmem;
    }
  return ajouter_(variable_imposee_,secmem);
}
 
DoubleTab& Source_PDF_base::ajouter_(const DoubleTab& variable, DoubleTab& resu) const
{
  const int i_traitement_special = 0 ;
  ajouter_(variable, resu, i_traitement_special) ;
  return resu;
}
 
DoubleTab& Source_PDF_base::ajouter_(const DoubleTab& variable, DoubleTab& resu, const int i_traitement_special) const
{
  Cerr << "Source_PDF_base: Not implemented for current discretisation. Aborting..." << finl;
  abort();
  return resu;
}
 
void  Source_PDF_base::contribuer_a_avec(const DoubleTab& inco, Matrice_Morse& matrice) const
{
  Cerr << "Source_PDF_base: Not implemented for current discretisation. Aborting..." << finl;
  abort();
}
 
DoubleTab& Source_PDF_base::calculer(DoubleTab& resu) const
{
  resu = 0;
  const DoubleTab& variable=equation().inconnue().valeurs();
  return ajouter_(variable,resu);
}
 
DoubleVect& Source_PDF_base::compute_source_term_PDF(DoubleVect& bilan)
{
  int i_traitement_special = 0;
  const DoubleTab& variable=equation().inconnue().valeurs();
  DoubleTrav add_zero_term(variable);
  add_zero_term = 0.;
  return compute_source_term_PDF(i_traitement_special, add_zero_term, bilan);
}
 
DoubleVect& Source_PDF_base::compute_source_term_PDF(int i_traitement_special, DoubleTab& add_term, DoubleVect& bilan)
{
  assert(Objet_U::dimension <= 3);
  // double temps_corrige = temps_computation_pdf_ + dt_computation_pdf_;
  // Cerr<<"compute_source_term_PDF:: current time, temps computation_pdf and corrected time = "<<equation().probleme().schema_temps().temps_courant() <<" "<<temps_computation_pdf_<<" "<<temps_corrige<<finl;
  DoubleTab& variable=equation().inconnue().valeurs();
  // If the current time has not been updated yet, use the future time for the unknown
  if (equation().probleme().schema_temps().temps_courant() == temps_computation_pdf_) variable=equation().inconnue().futur();
  int nb_som= variable.dimension(0);
  int nb_comp = variable.dimension(1);
  assert(nb_comp == modele_lu_.dim_variable_);
 
  DoubleTrav resu(variable);
  calculer(resu, i_traitement_special);
  // Make sure we have the same dt as for sec_mem_pdf (calculer_pdf)
  double dt_courant = calcul_dt_pdf();
  resu *= dt_courant/dt_computation_pdf_;
  DoubleTrav flag_cl(resu);
  filtre_CLD(flag_cl);
 
  source_term_PDF.resize(sec_mem_pdf.dimension(0), sec_mem_pdf.dimension(1));
  source_term_PDF = 0.;
 
  DoubleVect source_term[3];
  bilan.resize(equation().inconnue().nb_comp());
  bilan = 0.0;
  Cerr<<"(IBM) Source_PDF_base: Bilan terme PDF = ";
  for (int i=0; i<nb_comp; i++)
    {
      source_term[i] = champ_nodal_->valeurs();
      source_term[i] = 0.;
      for (int j=0; j<nb_som; j++)
        {
          double  filter = (std::fabs(resu(j,i)) > 0?1.:0.);
          filter *= flag_cl(j,i);
          source_term[i](j) = (resu(j,i) - sec_mem_pdf(j,i) + add_term(j,i) )*filter;
          source_term_PDF(j,i) = source_term[i](j);
          // if(  (abs(resu(j,i)) > 1.0) )
          //   {
          //     Cerr<<"i,j = "<<i<<" "<<j<<" ; resu = "<<resu(j,i)<<" ; sec_mem_pdf = "<<sec_mem_pdf(j,i)<<" ; bilan = "<<source_term[i](j)<<finl;
          //     Cerr<<"added term (secmem_conv) = "<<add_term(j,i)<<finl;
          //   }
        }
      bilan(i) = mp_somme_vect(source_term[i]);
      Cerr<< bilan(i)<<" ";
    }
  source_term_PDF.echange_espace_virtuel();
  Cerr<<finl;
  return bilan;
}
 
void Source_PDF_base::volume_source_term_PDF(DoubleTab& volume_term_PDF)
{
}
 
DoubleTab& Source_PDF_base::calculer(DoubleTab& resu, const int i_traitement_special) const
{
  resu = 0;
  const DoubleTab& variable=equation().inconnue().valeurs();
  return ajouter_(variable,resu, i_traitement_special);
}
 
void Source_PDF_base::calculer_variable_imposee_elem_fluid()
{
  Cerr << "Source_PDF_base: Not implemented for current discretisation. Aborting..." << finl;
  abort();
}
 
void Source_PDF_base::calculer_variable_imposee_mean_grad()
{
  Cerr << "Source_PDF_base: Not implemented for current discretisation. Aborting..." << finl;
  abort();
}
 
void Source_PDF_base::calculer_variable_imposee_hybrid()
{
  Cerr << "Source_PDF_base: Not implemented for current discretisation. Aborting..." << finl;
  abort();
}
 
void Source_PDF_base::calculer_vitesse_imposee_power_law_tbl()
{
  Cerr << "Source_PDF_base: Not implemented for current discretisation. Aborting..." << finl;
  abort();
}
 
void Source_PDF_base::calculer_vitesse_imposee_power_law_tbl_u_star()
{
  Cerr << "Source_PDF_base: Not implemented for current discretisation. Aborting..." << finl;
  abort();
}
 
void Source_PDF_base::calculer_temperature_imposee_wall_law()
{
  Cerr << "Source_PDF_base: Not implemented for current discretisation. Aborting..." << finl;
  abort();
}
 
void Source_PDF_base::calculer_variable_imposee()
{
  if (interpolation_bool_)
    {
      if (interpolation_lue_->que_suis_je() == "Interpolation_IBM_element_fluide")
        calculer_variable_imposee_elem_fluid();
      else if (interpolation_lue_->que_suis_je() == "Interpolation_IBM_gradient_moyen")
        calculer_variable_imposee_mean_grad();
      else if (interpolation_lue_->que_suis_je() == "Interpolation_IBM_hybride")
        calculer_variable_imposee_hybrid();
      else if (interpolation_lue_->que_suis_je() == "Interpolation_IBM_power_law_tbl")
        calculer_vitesse_imposee_power_law_tbl();
      else if (interpolation_lue_->que_suis_je() == "Interpolation_IBM_power_law_tbl_u_star")
        calculer_vitesse_imposee_power_law_tbl_u_star();
      else if (interpolation_lue_->que_suis_je() == "Interpolation_IBM_thermal_wall_law")
        calculer_temperature_imposee_wall_law();
    }
}
 
DoubleTab& Source_PDF_base::calculer_pdf(DoubleTab& resu) const
{
  resu = 0;
  ajouter_(variable_imposee_,resu);
 
  double my_dt = (ref_cast_non_const(Source_PDF_base, *this)).calcul_dt_pdf();
  (ref_cast_non_const(Source_PDF_base, *this)).set_dt_computation_pdf(my_dt);
  double my_temps = equation().probleme().schema_temps().temps_courant();
  (ref_cast_non_const(Source_PDF_base, *this)).set_temps_computation_pdf(my_temps);
 
  return resu;
}
 
double Source_PDF_base::calcul_dt_pdf()
{
  double dt_ref = equation().probleme().schema_temps().pas_de_temps();
  double dt_min = equation().probleme().schema_temps().pas_temps_min();
  double my_dt = std::max(dt_ref,dt_min);
  if (equation().probleme().schema_temps().temps_courant()==0) my_dt = 1.;
  return my_dt;
}
 
void Source_PDF_base::mettre_a_jour(double temps)
{
  //la_source->mettre_a_jour(temps);
}
 
void Source_PDF_base::correct_variable(const DoubleTab& coeff_node, DoubleTab& variable) const
{
  int nb_data_vit = variable.size();
  int nb_node_2 = coeff_node.dimension(0) ;
  int ncomp = coeff_node.dimension(1) ;
  if (nb_data_vit != nb_node_2 * ncomp)
    {
      Cerr<<"Source_PDF_base: numbers of nodes are different for variable and coeff_node."<<finl;
      abort();
    }
  const DoubleTab& vit_dir=variable_imposee_;
  for(int i=0; i<nb_node_2; i++)
    {
      for (int j=0; j<ncomp; j++)
        {
          if (coeff_node(i,j)!= 1.0)
            {
              variable(i,j)=vit_dir(i,j);
            }
          // Cerr<<"Source_PDF_base: coeff_node variable("<<i<<","<<j<<"): "<<coeff_node(i,j)<<" "<<variable(i,j)<<finl;
        }
    }
}
 
int Source_PDF_base::impr(Sortie& os) const
{
  Cerr << "Source_PDF_base: Not implemented for current discretisation. Aborting..." << finl;
  abort();
  return 0;
}
 
void Source_PDF_base::ouvrir_fichier(SFichier& os, const Nom& type, const int flag=1) const
{
  // null flag: do not open the file
  if (flag==0)
    return ;
 
  const Probleme_base& pb=equation().probleme();
  const Schema_Temps_base& sch=pb.schema_temps();
  const int precision=sch.precision_impr();
  Nom nomfichier(out_);
  if (type!="") nomfichier+=(Nom)"_"+type;
  nomfichier+=".out";
 
  // Create the file at the first output with the header
  if (sch.nb_impr()==1 && !pb.reprise_effectuee())
    {
      os.ouvrir(nomfichier);
      SFichier& fic=os;
      Nom espace="\t\t";
      fic << (Nom)"# Printing of the source term "+que_suis_je()+" of the equation "+equation().que_suis_je()+" of the problem "+equation().probleme().le_nom() << finl;
      fic << "# " << description() << finl;
      if(modele_lu_.dim_variable_ == Objet_U::dimension)
        {
          // vector
          assert(Objet_U::dimension==3);
          fic << "# Time" << espace << "Fx" << espace << "Fy" << espace << "Fz";
        }
      else if(modele_lu_.dim_variable_ == 1)
        {
          // scalar
          fic << "# Time" << espace << "Sum";
        }
      else
        {
          Cerr << "Source_PDF_base: for scalar or vector only; dim = " << modele_lu_.dim_variable_ << finl;
          Process::exit();
        }
      fic << finl;
    }
  // Sinon on l'ouvre
  else
    os.ouvrir(nomfichier,ios::app);
 
  os.precision(precision);
  os.setf(ios::scientific);
}
 
void Source_PDF_base::updateChampRho()
{
  if (equation().probleme().que_suis_je() == "Pb_Melange")
    champ_rho_->affecter(equation().probleme().get_champ("masse_volumique_melange"));
  else
    champ_rho_->affecter(equation().probleme().get_champ("masse_volumique"));
}
 
void Source_PDF_base::correct_incr_pressure(const DoubleTab& coeff_node, DoubleTab& correction_en_pression) const
{
  Cerr << "Source_PDF_NS_base: Not implemented for current discretisation. Aborting..." << finl;
  abort();
}
 
void Source_PDF_base::correct_pressure(const DoubleTab& coeff_node, DoubleTab& pression, const DoubleTab& correction_en_pression) const
{
  Cerr << "Source_PDF_NS_base: Not implemented for current discretisation. Aborting..." << finl;
  abort();
}
 
void Source_PDF_base::creer_champ(const Motcle& motlu)
{
  const Probleme_base& pb = equation().probleme();
  const DoubleTab& variable=equation().inconnue().valeurs();
  int nb_comp = variable.dimension(1);
  Motcle directive("temperature");
  if (equation().inconnue().is_vectorial()) directive="vitesse";
  const Domaine_dis_base& le_dom_dis = equation().domaine_dis();
 
  double temps=0.;
  if (motlu=="source_term_pdf" && !champ_source_term_PDF_)
    {
      if (nb_comp == 1)
        {
          Noms noms(1);
          noms[0]="source_term_PDF";
          Noms unites(1);
          unites[0] = "-";
          pb.discretisation().discretiser_champ(directive,le_dom_dis,scalaire,noms,unites,nb_comp,temps,champ_source_term_PDF_);
        }
      else if (nb_comp == Objet_U::dimension)
        {
          Noms noms(nb_comp);
          noms[0]="source_term_PDF";
          Noms unites(nb_comp);
          unites[0] = "-";
          pb.discretisation().discretiser_champ(directive,le_dom_dis,vectoriel,noms,unites,nb_comp,temps,champ_source_term_PDF_);
        }
      champs_compris_.ajoute_champ(champ_source_term_PDF_);
    }
 
  if (motlu=="barycentre_IBM" && !champ_barycentre_IBM_)
    {
      nb_comp = Objet_U::dimension;
      Noms nomsb(nb_comp);
      nomsb[0]="barycentre_IBM";
      Noms unitesb(nb_comp);
      unitesb[0] = "m";
      pb.discretisation().discretiser_champ("champ_elem",le_dom_dis,vectoriel,nomsb,unitesb,nb_comp,0., champ_barycentre_IBM_);
      champs_compris_.ajoute_champ(champ_barycentre_IBM_);
    }
 
  if (motlu=="normal_IBM" && !champ_normal_IBM_)
    {
      nb_comp = Objet_U::dimension;
      Noms nomsn(nb_comp);
      nomsn[0]="normal_IBM";
      Noms unitesn(nb_comp);
      unitesn[0] = "m";
      pb.discretisation().discretiser_champ("champ_elem",le_dom_dis,vectoriel,nomsn,unitesn,nb_comp,0., champ_normal_IBM_);
      champs_compris_.ajoute_champ(champ_normal_IBM_);
    }
 
  if (motlu=="aire_IBM" && !champ_aire_IBM_)
    {
      pb.discretisation().discretiser_champ("champ_elem",le_dom_dis,"aire_IBM","m^2",1,0., champ_aire_IBM_);
      champs_compris_.ajoute_champ(champ_aire_IBM_);
    }
 
  if (motlu=="velocity_shape_IBM" && !champ_vitesse_shape_IBM_)
    {
      nb_comp = Objet_U::dimension;
      Noms nomsv(nb_comp);
      for (int i=0; i< nb_comp; i++) nomsv[i]="velocity_shape_IBM";
      Noms unitesv(nb_comp);
      for (int i=0; i< nb_comp; i++) unitesv[i] = "m/s";
      pb.discretisation().discretiser_champ("champ_elem",le_dom_dis,vectoriel,nomsv,unitesv,nb_comp,0., champ_vitesse_shape_IBM_);
      champs_compris_.ajoute_champ(champ_vitesse_shape_IBM_);
    }
 
  if (motlu=="pseudo_level_set_IBM" && !champ_pseudo_level_set_IBM_)
    {
      nb_comp = 1;
      Noms nompl(nb_comp);
      for (int i=0; i< nb_comp; i++) nompl[i]="pseudo_level_set_IBM";
      Noms unitepl(nb_comp);
      for (int i=0; i< nb_comp; i++) unitepl[i] = "m";
      Motcle directiveps("temperature");
      pb.discretisation().discretiser_champ(directive,le_dom_dis,scalaire,nompl,unitepl,nb_comp,0., champ_pseudo_level_set_IBM_);
      champs_compris_.ajoute_champ(champ_pseudo_level_set_IBM_);
    }
}
 
void Source_PDF_base::get_noms_champs_postraitables(Noms& nom,Option opt) const
{
  Noms noms_compris = champs_compris_.liste_noms_compris();
  noms_compris.add("source_term_pdf");
  noms_compris.add("barycentre_IBM");
  noms_compris.add("normal_IBM");
  noms_compris.add("aire_IBM");
  noms_compris.add("velocity_shape_IBM");
  noms_compris.add("pseudo_level_set_IBM");
  if (opt==DESCRIPTION)
    Cerr<<" Source_PDF_base : "<< noms_compris <<finl;
  else
    nom.add(noms_compris);
}
 
bool Source_PDF_base::has_champ(const Motcle& nom, OBS_PTR(Champ_base) &ref_champ) const
{
  if (nom == "source_term_pdf" && champ_source_term_PDF_)
    {
      ref_champ = get_champ(nom);
      return true;
    }
  else if (nom == "barycentre_IBM" && champ_barycentre_IBM_)
    {
      ref_champ = get_champ(nom);
      return true;
    }
  else if (nom == "normal_IBM" && champ_normal_IBM_)
    {
      ref_champ = get_champ(nom);
      return true;
    }
  else if (nom == "aire_IBM" && champ_aire_IBM_)
    {
      ref_champ = get_champ(nom);
      return true;
    }
  else if (nom == "velocity_shape_IBM" && champ_vitesse_shape_IBM_)
    {
      ref_champ = get_champ(nom);
      return true;
    }
  else if (nom == "pseudo_level_set_IBM" && champ_pseudo_level_set_IBM_)
    {
      ref_champ = get_champ(nom);
      return true;
    }
  else
    return false; /* nothing found */
}
 
bool Source_PDF_base::has_champ(const Motcle& nom) const
{
  if (nom == "source_term_pdf" && champ_source_term_PDF_)
    return true;
  else if (nom == "barycentre_IBM" && champ_barycentre_IBM_)
    return true;
  else if (nom == "normal_IBM" && champ_normal_IBM_)
    return true;
  else if (nom == "aire_IBM" && champ_aire_IBM_)
    return true;
  else if (nom == "velocity_shape_IBM" && champ_vitesse_shape_IBM_)
    return true;
  else if (nom == "pseudo_level_set_IBM" && champ_pseudo_level_set_IBM_)
    return true;
  else
    return false; /* nothing found */
}
 
const Champ_base& Source_PDF_base::get_champ(const Motcle& nom) const
{
  int dim_esp =  Objet_U::dimension;
 
  if (nom=="source_term_pdf")
    {
      if (!champ_source_term_PDF_)
        throw std::runtime_error(std::string("Field ") + nom.getString() + std::string(" not found !"));
 
      // Initialization to 0 of champ_source_term_PDF_
      DoubleTab& valeurs = champ_source_term_PDF_->valeurs();
      valeurs=0.;
      if (source_term_PDF.size_array()>0)
        {
          int nb_d0 = source_term_PDF.dimension(0);
          int nb_d1 = source_term_PDF.dimension(1);
          for (int num0=0; num0<nb_d0; num0++)
            for (int num1=0; num1<nb_d1; num1++)
              valeurs(num0,num1)=source_term_PDF(num0,num1);
        }
      valeurs.echange_espace_virtuel();
      champ_source_term_PDF_ ->mettre_a_jour(equation().probleme().schema_temps().temps_courant());
      return champs_compris_.get_champ(nom);
    }
  else if (nom=="barycentre_IBM")
    {
      if (!champ_barycentre_IBM_)
        throw std::runtime_error(std::string("Field ") + nom.getString() + std::string(" not found !"));
 
      // Initialization to 0 of champ_barycentre_IBM_
      DoubleTab& valeurs = champ_barycentre_IBM_->valeurs();
      valeurs=0.;
      if (champ_barycentre_)
        {
          const DoubleTab& barycentre = champ_barycentre_->valeurs();
          int nb_d0 = barycentre.dimension(0);
          int nb_d1 = barycentre.dimension(1);
          for (int num0=0; num0<nb_d0; num0++)
            for (int num1=0; num1<nb_d1; num1++)
              valeurs(num0,num1)=barycentre(num0,num1);
        }
      valeurs.echange_espace_virtuel();
      champ_barycentre_IBM_->mettre_a_jour(equation().probleme().schema_temps().temps_courant());
      return champs_compris_.get_champ(nom);
    }
  else if (nom=="normal_IBM")
    {
      if (!champ_normal_IBM_)
        throw std::runtime_error(std::string("Field ") + nom.getString() + std::string(" not found !"));
 
      // Initialization to 0 of champ_normal_IBM_
      DoubleTab& valeurs = champ_normal_IBM_->valeurs();
      valeurs=0.;
      if (champ_rotation_)
        {
          const DoubleTab& rotation = champ_rotation_->valeurs();
          int nb_d0 = rotation.dimension(0);
          for (int num0=0; num0<nb_d0; num0++)
            for (int num1=0; num1<dim_esp; num1++)
              valeurs(num0,num1)=rotation(num0,3*num1+(dim_esp-1));
        }
      valeurs.echange_espace_virtuel();
      champ_normal_IBM_->mettre_a_jour(equation().probleme().schema_temps().temps_courant());
      return champs_compris_.get_champ(nom);
    }
  else if (nom=="aire_IBM")
    {
      if (!champ_aire_IBM_)
        throw std::runtime_error(std::string("Field ") + nom.getString() + std::string(" not found !"));
 
      // Initialization to 0 of champ_aire_IBM_
      DoubleTab& valeurs = champ_aire_IBM_->valeurs();
      valeurs=0.;
      if (champ_aire_)
        {
          const DoubleTab& aire = champ_aire_->valeurs();
          int nb_d0 = aire.dimension(0);
          int nb_d1 = aire.dimension(1);
          for (int num0=0; num0<nb_d0; num0++)
            for (int num1=0; num1<nb_d1; num1++)
              valeurs(num0,num1)=aire(num0,num1);
        }
      valeurs.echange_espace_virtuel();
      champ_aire_IBM_->mettre_a_jour(equation().probleme().schema_temps().temps_courant());
      return champs_compris_.get_champ(nom);
    }
  else if (nom=="velocity_shape_IBM")
    {
      if (!champ_vitesse_shape_IBM_)
        throw std::runtime_error(std::string("Field ") + nom.getString() + std::string(" not found !"));
 
      // Initialization to 0 of champ_vitesse_shape_IBM_
      DoubleTab& valeurs = champ_vitesse_shape_IBM_->valeurs();
      valeurs=0.;
      if (modele_lu_.vitesse_shape_IBM_)
        {
          const DoubleTab& vitesse = modele_lu_.get_vitesse_shape_IBM();
          int nb_d0 = vitesse.dimension(0);
          int nb_d1 = vitesse.dimension(1);
          for (int num0=0; num0<nb_d0; num0++)
            for (int num1=0; num1<nb_d1; num1++)
              valeurs(num0,num1)=vitesse(num0,num1);
        }
      valeurs.echange_espace_virtuel();
      champ_vitesse_shape_IBM_->mettre_a_jour(equation().probleme().schema_temps().temps_courant());
      return champs_compris_.get_champ(nom);
    }
  else if (nom=="pseudo_level_set_IBM")
    {
      if (!champ_pseudo_level_set_IBM_)
        throw std::runtime_error(std::string("Field ") + nom.getString() + std::string(" not found !"));
 
      // Initialization to 0 of champ_pseudo_level_set_IBM_
      DoubleTab& valeurs = champ_pseudo_level_set_IBM_->valeurs();
      valeurs=0.;
      if (interpolation_bool_)
        {
          if (interpolation_lue_->pseudo_level_set_)
            {
              Interpolation_IBM_base& my_interp = ref_cast_non_const(Interpolation_IBM_base, getInterpolationLu());
              const DoubleTab& pseudo_level_set = my_interp.get_pseudo_level_set();
              int nb_d0 = pseudo_level_set.dimension(0);
              for (int num0=0; num0<nb_d0; num0++) valeurs(num0)=pseudo_level_set(num0);
            }
        }
      valeurs.echange_espace_virtuel();
      champ_pseudo_level_set_IBM_->mettre_a_jour(equation().probleme().schema_temps().temps_courant());
      return champs_compris_.get_champ(nom);
    }
  else
    return champs_compris_.get_champ(nom);
}