Champ_front_recyclage.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 <Champ_front_recyclage.h>
#include <Interprete.h>
#include <Probleme_base.h>
#include <TRUSTTrav.h>
#include <Domaine_VF.h>
#include <Equation_base.h>
#include <communications.h>
#include <Octree.h>
#include <TRUSTList.h>
#include <EFichier.h>
#include <LecFicDiffuse.h>
#include <EcrFicCollecte.h>
#include <Champ_front_calc.h>
#include <string> // For AIX
#include <Param.h>
 
Implemente_instanciable_sans_constructeur(Champ_front_recyclage,"Champ_front_recyclage",Ch_front_var_instationnaire_dep);
// XD champ_front_recyclage front_field_base champ_front_recyclage INHERITS_BRACE This keyword is used on a boundary to
// XD_CONT get a field from another boundary. NL2 It is to use, in a general way, on a boundary of a local_pb problem, a
// XD_CONT field calculated from a linear combination of an imposed field g(x,y,z,t) with an instantaneous f(x,y,z,t)
// XD_CONT and a spatial mean field <f>(t) or a temporal mean field <f>(x,y,z) extracted from a plane of a problem named
// XD_CONT pb (pb may be local_pb itself): NL2 For each component i, the field F applied on the boundary will be: NL2
// XD_CONT F_i(x,y,z,t) = alpha_i*g_i(x,y,z,t) + xsi_i*[f_i(x,y,z,t)- beta_i*<fi>]
// XD attr pb_champ_evaluateur pb_champ_evaluateur pb_champ_evaluateur REQ not_set
// XD attr distance_plan listf distance_plan OPT Vector which gives the distance between the boundary and the plane from
// XD_CONT where the field F will be extracted. By default, the vector is zero, that should imply the two domains have
// XD_CONT coincident boundaries.
// XD attr ampli_moyenne_imposee list ampli_moyenne_imposee OPT 2|3 alpha(0) alpha(1) [alpha(2)]: alpha_i coefficients
// XD_CONT (by default =1)
// XD attr ampli_moyenne_recyclee list ampli_moyenne_recyclee OPT 2|3 beta(0) beta(1) [beta(2)]}: beta_i coefficients
// XD_CONT (by default =1)
// XD attr ampli_fluctuation list ampli_fluctuation OPT 2|3 gamma(0) gamma(1) [gamma(2)]}: gamma_i coefficients (by
// XD_CONT default =1)
// XD attr direction_anisotrope entier(into=[1,2,3]) direction_anisotrope OPT If an integer is given for direction (X:1,
// XD_CONT Y:2, Z:3, by default, direction is negative), the imposed field g will be 0 for the 2 other directions.
// XD attr moyenne_imposee moyenne_imposee_deriv moyenne_imposee OPT Value of the imposed g field.
// XD attr moyenne_recyclee chaine moyenne_recyclee OPT Method used to perform a spatial or a temporal averaging of f
// XD_CONT field to specify <f>. <f> can be the surface mean of f on the plane (surface option, see below) or it can be
// XD_CONT read from several files (for example generated by the chmoy_faceperio option of the Traitement_particulier
// XD_CONT keyword to obtain a temporal mean field). The option methode_recyc can be: surfacique, Surface mean for <f>
// XD_CONT from f values on the plane ; Or one of the following methode_moy options applied to read a temporal mean
// XD_CONT field <f>(x,y,z): interpolation, connexion_approchee or connexion_exacte
// XD attr fichier chaine fichier OPT not_set
 
// XD pb_champ_evaluateur objet_u nul NO_BRACE specifies problem name, the field name beloging to the problem and number
// XD_CONT of field components.
// XD attr pb ref_Pb_base pb REQ name of the problem where the source fields will be searched.
// XD attr champ chaine champ REQ name of the field
// XD attr ncomp entier ncomp REQ number of components
 
Champ_front_recyclage::Champ_front_recyclage()
{
  int dim = Objet_U::dimension;
  delt_dist.resize(dim,RESIZE_OPTIONS::NOCOPY_NOINIT);
  delt_dist=0;
  ampli_fluct_.resize(1,RESIZE_OPTIONS::NOCOPY_NOINIT);
  ampli_fluct_ = 1.;
  ampli_moy_imposee_.resize(1,RESIZE_OPTIONS::NOCOPY_NOINIT);
  ampli_moy_imposee_ = 1.;
  ampli_moy_recyclee_.resize(1,RESIZE_OPTIONS::NOCOPY_NOINIT);
  ampli_moy_recyclee_ = 1.;
 
  methode_moy_impos_ = -1;
  methode_moy_recycl_ = -1;
  ndir = -4;
  diametre = -1.;
  u_tau = -1.;
  visco_cin = -1.;
  fich_recycl_ = "??";
  fich_impos_ = "??";
}
 
Sortie& Champ_front_recyclage::printOn(Sortie& os) const
{
  return os;
}
 
//-pb_champ_evaluateur is mandatory to specify the name of the problem
// and the name of the unknown that will be used to initialize the reference
// to the evaluator field. The number of components of this field must be specified.
 
Entree& Champ_front_recyclage::readOn(Entree& is)
{
  Param param(que_suis_je());
  set_param(param);
  param.lire_avec_accolades_depuis(is);
  return is;
}
 
void Champ_front_recyclage::set_param(Param& param) const
{
  param.ajouter_non_std("pb_champ_evaluateur",(this),Param::REQUIRED);
  param.ajouter_non_std("distance_plan",(this));
  param.ajouter_non_std("ampli_fluctuation",(this));
  param.ajouter_non_std("ampli_moyenne_imposee",(this));
  param.ajouter_non_std("ampli_moyenne_recyclee",(this));
  param.ajouter_non_std("moyenne_imposee",(this));
  param.ajouter_non_std("moyenne_recyclee",(this));
  param.ajouter("direction_anisotrope",&ndir);
}
 
int Champ_front_recyclage::lire_motcle_non_standard(const Motcle& mot, Entree& is)
{
  int retval = 1;
  if (mot=="pb_champ_evaluateur")
    {
      int nbcomp;
      is >> nom_pb1 >> nom_inco1 >> nbcomp;
      ampli_fluct_.resize(nbcomp,RESIZE_OPTIONS::NOCOPY_NOINIT);
      ampli_moy_imposee_.resize(nbcomp,RESIZE_OPTIONS::NOCOPY_NOINIT);
      ampli_moy_recyclee_.resize(nbcomp,RESIZE_OPTIONS::NOCOPY_NOINIT);
      ampli_fluct_ = 1.;
      ampli_moy_imposee_ = 1.;
      ampli_moy_recyclee_ = 1.;
      fixer_nb_comp(nbcomp);
    }
  else if (mot=="distance_plan")
    {
      for (int i=0; i<dimension; i++) is >> delt_dist(i) ;
    }
  else if (mot=="ampli_fluctuation")
    {
      int nbcomp;
      is >> nbcomp;
      ampli_fluct_.resize(nbcomp);
      for (int i=0; i<nbcomp; i++) is >> ampli_fluct_(i) ;
      fixer_nb_comp(nbcomp);
    }
  else if (mot=="ampli_moyenne_imposee")
    {
      int nbcomp;
      is >> nbcomp;
      ampli_moy_imposee_.resize(nbcomp);
      for (int i=0; i<nbcomp; i++) is >> ampli_moy_imposee_(i) ;
      fixer_nb_comp(nbcomp);
    }
  else if (mot=="ampli_moyenne_recyclee")
    {
      int nbcomp;
      is >> nbcomp;
      ampli_moy_recyclee_.resize(nbcomp);
      for (int i=0; i<nbcomp; i++) is >> ampli_moy_recyclee_(i) ;
      fixer_nb_comp(nbcomp);
    }
  else if (mot=="moyenne_imposee")
    {
      methode_moy_impos_ = lire_info_moyenne_imposee(is);
    }
  else if (mot=="moyenne_recyclee")
    {
      methode_moy_recycl_ = lire_info_moyenne_recyclee(is);
    }
  else retval = -1;
 
  return retval;
}
 
//Reading the information necessary to evaluate the imposed mean
//-profil                : to impose an analytical profile
//-interpolation       : reads from a file and builds a mean field
//                       by interpolating the data read.
//                       The mean is built for a preferred direction
//                         (direction_anisotrope) and is 0 for the other directions
//-connexion_approchee : reads from a file and retains the value of the
//                       variable read by connection with the closest point
//                       to the considered boundary face
//-connexion_exacte    : reads a geometry file containing the coordinates of points
//                         located in the evaluation plane and reads from a separate
//                       file the mean values. The mean values read
//                         are stored when the exact correspondence between
//                         facing points is verified.
//- logarithmique      : builds the mean using a wall law (logarithmic)
//
 
// XD moyenne_imposee_deriv objet_u moyenne_imposee_deriv INHERITS_BRACE not_set
 
// XD moyenne_imposee_profil moyenne_imposee_deriv profil NO_BRACE To specify analytic profile for the imposed g field.
// XD attr profile listchaine profile REQ specifies the analytic profile: 2|3 valx(x,y,z,t) valy(x,y,z,t)
// XD_CONT [valz(x,y,z,t)]
 
// XD moyenne_imposee_connexion_exacte moyenne_imposee_deriv connexion_exacte NO_BRACE To read the imposed field from
// XD_CONT two files.
// XD attr fichier chaine(into=["fichier"]) fichier REQ not_set
// XD attr file1 chaine file1 REQ first file, contains the points coordinates (which should be the same as the
// XD_CONT coordinates of the boundary faces). The format of this file is:NL2 N NL2 1 x(1) y(1) [z(1)] NL2 2 x(2) y(2)
// XD_CONT [z(2)] NL2 ... NL2 N x(N) y(N) [z(N)]
// XD attr file2 chaine file2 OPT second file, contains the mean values. The format of this file is: NL2 N NL2 1 valx(1)
// XD_CONT valy(1) [valz(1)] NL2 2 valx(2) valy(2) [valz(2)] NL2 ... NL2 N valx(N) valy(N) [valz(N)]
 
// XD moyenne_imposee_connexion_approchee moyenne_imposee_deriv connexion_approchee NO_BRACE To read the imposed field
// XD_CONT from a file where positions and values are given (it is not necessary that the coordinates of points match
// XD_CONT the coordinates of the boundary faces, indeed, the nearest point of each face of the boundary will be used).
// XD attr fichier chaine(into=["fichier"]) fichier REQ not_set
// XD attr file1 chaine file1 REQ filename. The format of the file is: NL2 N NL2 x(1) y(1) [z(1)] valx(1) valy(1)
// XD_CONT [valz(1)] NL2 x(2) y(2) [z(2)] valx(2) valy(2) [valz(2)] NL2 ... NL2 x(N) y(N) [z(N)] valx(N) valy(N)
// XD_CONT [valz(N)]
 
// XD moyenne_imposee_interpolation moyenne_imposee_deriv interpolation NO_BRACE To create an imposed field built by
// XD_CONT interpolation of values read from a file. The imposed field is applied on the direction given by the keyword
// XD_CONT direction_anisotrope (the field is zero for the other directions).
// XD attr fichier chaine(into=["fichier"]) fichier REQ The format of the file is: NL2 pos(1) val(1) NL2 pos(2) val(2)
// XD_CONT NL2 ... NL2 pos(N) val(N) NL2 If direction given by direction\_anisotrope is 1 (or 2 or 3), then pos will be
// XD_CONT X (or Y or Z) coordinate and val will be X value (or Y value, or Z value) of the imposed field.
// XD attr file1 chaine file1 REQ name of geom_face_perio
 
// XD moyenne_imposee_logarithmique moyenne_imposee_deriv logarithmique NO_BRACE To specify the imposed field (in this
// XD_CONT case, velocity) by an analytical logarithmic law of the wall: NL2 g(x,y,z) = u_tau * (
// XD_CONT log(0.5*diametre*u_tau/visco_cin)/Kappa + 5.1 ) NL2 with g(x,y,z)=u(x,y,z) if direction is set to 1,
// XD_CONT g=v(x,y,z) if direction is set to 2 and g=w(w,y,z) if it is set to 3
// XD attr diametre chaine(into=["diametre"]) diametre REQ not_set
// XD attr val floattant val REQ diameter
// XD attr u_tau chaine(into=["u_tau"]) u_tau REQ not_set
// XD attr val_u_tau floattant val_u_taul REQ value of u_tau
// XD attr visco_cin chaine(into=["visco_cin"]) visco_cin REQ not_set
// XD attr val_visco_cin floattant val_visco_cin REQ value of visco_cin
// XD attr direction chaine(into=["direction"]) direction REQ not_set
// XD attr val_direction entier val_direction REQ direction
 
int Champ_front_recyclage::lire_info_moyenne_imposee(Entree& is)
{
  int methode = -1;;
  Motcle choix,motlu;
  Motcles les_mots(5);
  {
    les_mots[0] = "profil";
    les_mots[1] = "interpolation";
    les_mots[2] = "connexion_approchee";
    les_mots[3] = "connexion_exacte";
    les_mots[4] = "logarithmique";
  }
 
  is >> choix;
  int rang = les_mots.search(choix);
  switch(rang)
    {
 
    case 0 :
      {
        methode = 1;
        int nc;
        is >> nc;
        fcts_xyz.dimensionner(nc);
        for (int i=0; i<nc; i++)
          is >> fcts_xyz[i];
        break;
      }
    case 1 :
      {
        methode = 2;
        break;
      }
    case 2 :
      {
        methode = 3;
        break;
      }
    case 3 :
      {
        methode = 4;
        break;
      }
    case 4 :
      {
        int compt_mot_lu = 0;
        while (compt_mot_lu!=4)
          {
            is >> motlu;
            if (motlu=="diametre")
              is >> diametre;
            else if (motlu=="u_tau")
              is >> u_tau;
            else if (motlu=="visco_cin")
              is >> visco_cin;
            else if (motlu=="direction")
              is >> ndir;
            else
              {
                Cerr<<"The keyword "<<motlu<<" is not recognized."<<finl;
              }
            compt_mot_lu += 1;
          }
        if (diametre<=0.)
          {
            Cerr<<"Champ_front_recyclage : a positive diameter must be specified when a logarithmique injected averaged is used."<<finl;
            exit();
          }
        if (u_tau<=0.)
          {
            Cerr<<"Champ_front_recyclage : a positive shear velocity must be specified when a logarithmique injected averaged is used."<<finl;
            exit();
          }
        if (visco_cin<=0.)
          {
            Cerr<<"Champ_front_recyclage : a positive kinematic viscosity must be specified when a logarithmique injected averaged is used."<<finl;
            exit();
          }
        if (ndir<-3)
          {
            Cerr<<"Champ_front_recyclage : the direction must be specified when a logarithmique injected averaged is used."<<finl;
            exit();
          }
        methode = 5;
        break;
      }
 
    default :
      {
        Cerr<<"Champ_front_recyclage::lire_info_moyenne_imposee : keyword "<<choix<<" is not recognized."<<finl;
        Cerr<<"The recognized keywords are :"<<les_mots<<finl;
        exit();
      }
    }
 
  if ((methode == 2) || (methode == 3))
    {
      Nom nom_fich;
      is >> motlu;
      if (motlu!="fichier")
        {
          Cerr<<"Champ_front_recyclage::lire_info_moyenne_imposee : at least one file must be provided for this method."<<finl;
          exit();
        }
      else
        is >> nom_fich;
      fich_impos_ = nom_fich;
    }
 
  if (methode == 4)
    {
      Nom nom_fich;
      is >> motlu;
      if (motlu!="fichier")
        {
          Cerr<<"Champ_front_recyclage::lire_info_moyenne_imposee : at least one file must be provided for this method."<<finl;
          exit();
        }
      else
        is >> nom_fich;
      fich_maillage_ = nom_fich;
      is >> nom_fich;
      fich_impos_ = nom_fich;
    }
 
  if (methode==-1)
    {
      Cerr<<"Champ_front_recyclage::lire_info_moyenne_imposee : no method to make average evaluation has been given."<<finl;
      Cerr<<"Choose among the following keywords : les_mots"<<finl;
      exit();
    }
 
  return methode;
}
 
//Reading the information necessary to evaluate the recycled mean
//-surfacique : surface average of the recycled values
//             (the average is computed on bord2 where the values are retrieved)
//-interpolation connexion_approchee connexion_exacte (see lire_info_moyenne_imposee)
//
int Champ_front_recyclage::lire_info_moyenne_recyclee(Entree& is)
{
  int methode = -1;
  Motcle choix,motlu;
  Motcles les_mots(4);
  {
    les_mots[0] = "surfacique";
    les_mots[1] = "interpolation";
    les_mots[2] = "connexion_approchee";
    les_mots[3] = "connexion_exacte";
  }
  is >> choix;
  int rang = les_mots.search(choix);
 
  switch(rang)
    {
    case 0 :
      {
        methode = 1;
        break;
      }
    case 1 :
      {
        methode = 2;
        break;
      }
    case 2 :
      {
        methode = 3;
        break;
      }
    case 3 :
      {
        methode = 4;
        break;
      }
 
    default :
      {
        Cerr<<"Champ_front_recyclage::lire_info_moyenne_recyclee keyword "<<choix<<" is not recognized."<<finl;
        Cerr<<"The recognized keywords are :"<<les_mots<<finl;
        exit();
      }
    }
 
  if ((methode == 2) || (methode == 3) || (methode == 4))
    {
      Nom nom_fich;
      is >> motlu;
      if (motlu!="fichier")
        {
          Cerr<<"Champ_front_recyclage::lire_info_moyenne_recyclee : a file must be provided for this method."<<finl;
          exit();
        }
      else
        is >> nom_fich;
 
      fich_recycl_  = nom_fich;
      if (methode == 4)
        is >> nom_fich;
      fich_maillage_ = nom_fich;
    }
 
  if (methode==-1)
    {
      Cerr<<"Champ_front_recyclage::lire_info_moyenne_recyclee : no method to make average evaluation has been given."<<finl;
      Cerr<<"Choose among the following keywords : les_mots"<<finl;
      exit();
    }
 
  return methode;
}
 
//Initialization of the reference to the evaluator field (l_inconnue1)
void Champ_front_recyclage::associer_champ_evaluateur(const Nom& un_nom_pb1,const Motcle& un_nom_inco1)
{
  Objet_U& ob1 = Interprete::objet(un_nom_pb1);
  OBS_PTR(Probleme_base) pb1;
  OBS_PTR(Champ_base) rch1;
 
  if(sub_type(Probleme_base,ob1))
    {
      pb1 = ref_cast(Probleme_base,ob1);
    }
  else
    {
      Cerr <<"No problem named "<< un_nom_pb1 <<" has been found."<< finl;
      exit();
    }
  rch1 = pb1->get_champ(un_nom_inco1);
 
  if(sub_type(Champ_Inc_base,rch1.valeur()))
    l_inconnue1 = ref_cast(Champ_Inc_base,rch1.valeur()) ;
  else
    {
      Cerr<<pb1->le_nom()<<" has no unknown field named "<<un_nom_inco1<<finl;
      exit();
    }
}
 
//Evaluation of the location of points (coords)
//where l_inconnue1 must evaluate the values to be recycled
//coords = xv2(faces_bor2) + delt_dist
//
void Champ_front_recyclage::get_coord_faces(const Frontiere_dis_base& fr_vf,
                                            DoubleTab& coords,
                                            const DoubleVect& delt_dist)
{
  const Front_VF&   fr_vf2    = ref_cast(Front_VF,fr_vf);
  const Domaine_VF&    le_dom2  = ref_cast(Domaine_VF,fr_vf2.domaine_dis());
  const DoubleTab&  xv2       = le_dom2.xv();
  const int      nb_faces2 = fr_vf2.nb_faces();
  const int      ndeb2     = fr_vf2.num_premiere_face();
  const int      dim       = xv2.dimension(1);
 
  coords.resize(nb_faces2,dim);
  for (int i = 0; i < nb_faces2; i++)
    for (int j = 0; j < dim; j++)
      coords(i,j) = xv2(i+ndeb2,j) + delt_dist(j);
}
 
//-call to Champ_front_calc::initialiser(...)
//-dimensioning of moyenne_imposee_ and moyenne_recyclee_ by nb_faces_bord2 and nb_compo_
//-construction of the structures ensuring parallelism:
// inconnues1_coords_to_eval_ inconnues1_elems_ inconnues2_faces_
//-initialization of moyenne_imposee_ and moyenne_recyclee_
//
int Champ_front_recyclage::initialiser(double temps, const Champ_Inc_base& inco)
{
  if (!Ch_front_var_instationnaire_dep::initialiser(temps,inco))
    return 0;
 
  // 1: remote
  // 2: local
  const Front_VF& fr_vf2 = ref_cast(Front_VF,la_frontiere_dis.valeur());
  const Nom& nom_bord2 = fr_vf2.frontiere().le_nom();
  int nb_faces_bord2 = fr_vf2.nb_faces();
  associer_champ_evaluateur(nom_pb1,nom_inco1);
 
  moyenne_imposee_.resize(nb_faces_bord2,nb_compo_,RESIZE_OPTIONS::COPY_INIT);
  moyenne_recyclee_.resize(nb_faces_bord2,nb_compo_,RESIZE_OPTIONS::COPY_INIT);
 
  if (ampli_fluct_.size()!=nb_compo_)
    {
      Cerr<<"The number of fluctuation amplification factors is not coherent with the component number of the field "<<nom_inco1<<finl;
      Cerr<<"This field has "<<nb_compo_<<" components and there is "<<ampli_fluct_.size()<<" amplification factors defined."<<finl;
      Cerr<<"Please check that ampli_fluctuation is well specified in your data set."<<finl;
      exit();
    }
 
  if (ampli_moy_imposee_.size()!=nb_compo_)
    {
      Cerr<<"The number of the amplification factors for the fixed average is not coherent with the component number of the field "<<nom_inco1<<finl;
      Cerr<<"This field has "<<nb_compo_<<" components and there is "<<ampli_moy_imposee_.size()<<" amplification factors defined."<<finl;
      Cerr<<"Please check that ampli_moyenne_imposee is well specified in your data set."<<finl;
      exit();
    }
 
  if (ampli_moy_recyclee_.size()!=nb_compo_)
    {
      Cerr<<"The number of the amplification factors for the recycled average is not coherent with the component number of the field "<<nom_inco1<<finl;
      Cerr<<"This field has "<<nb_compo_<<" components and there is "<<ampli_moy_recyclee_.size()<<" amplification factors defined."<<finl;
      Cerr<<"Please check that ampli_moyenne_recyclee is well specified in your data set."<<finl;
      exit();
    }
 
  Cerr << "Initializing Champ_front_recyclage on the boundary: " << nom_bord2 << finl;
 
  // For each point in remote_coord, find the processor whose domain1 contains this point,
  //  and the index of the element:
 
  const int nprocs = nproc();
  const int moi = me();
 
  // Filling remote_coords[i]: coordinates of the face centres of boundary nom_bord2
  // (destination boundary) + delt_dist
  // of processor i (all procs hold the same remote_coords array)
  // Build a DoubleTabs containing the coordinates of the remote faces on domaine 1 for each processor
  DoubleTabs remote_coords(nprocs);
  /*
  if (nom_bord!="??") // Improve this test
  {
     const Front_VF& fr_vf1 = ref_cast(Front_VF,front_dis());
     get_coord_faces(fr_vf1,remote_coords[moi],delt_dist);
  }
  else */
  get_coord_faces(fr_vf2,remote_coords[moi],delt_dist);
 
  for (int pe = 0; pe < nprocs; pe++)
    envoyer_broadcast(remote_coords[pe],pe);
 
  // Size the 3 arrays:
  inconnues1_coords_to_eval_.dimensionner(nprocs);
  inconnues1_elems_.dimensionner(nprocs);
  inconnues2_faces_.dimensionner(nprocs);
 
  ArrsOfInt indexes_to_recv(nprocs);
 
  // tmp array for octree results:
  ArrOfInt elem_list;
 
 
  const Domaine& domaine1 = l_inconnue1->equation().domaine_dis().domaine();
  const int nb_elem_domaine1 = domaine1.nb_elem();
  const int dim = remote_coords[moi].dimension(1);
 
  const Domaine_dis_base& zdis = l_inconnue1->equation().domaine_dis();
  const Domaine_VF& zvf = ref_cast(Domaine_VF,zdis);
  const DoubleTab& xp = zvf.xp();
  DoubleVect remote_point(3);
  const OctreeRoot& octree = domaine1.construit_octree();
  // Loop on processes
  for (int pe = 0; pe < nprocs; pe++)
    {
      // Remote coordinates of the local faces:
      const DoubleTab& remote_coords_on_pe = remote_coords[pe];
      int nb_faces_on_pe = remote_coords_on_pe.dimension(0);
      elem_list.resize_tab(nb_faces_on_pe,RESIZE_OPTIONS::NOCOPY_NOINIT);
 
      ArrOfInt elem_list_par_pos;
      // Loop on local faces on the process pe:
      for (int ind_face=0; ind_face<nb_faces_on_pe; ind_face++)
        {
          for (int j=0; j<dimension; j++)
            remote_point(j) = remote_coords_on_pe(ind_face,j);
 
          // Fill elem_list_par_pos by elements surrounding the remote point:
          octree.rang_elems_sommet(elem_list_par_pos,remote_point(0),remote_point(1),remote_point(2));
          int size_elem_list = elem_list_par_pos.size_array();
 
          if (size_elem_list==0)
            elem_list[ind_face] = -1;
          else
            {
              DoubleTrav coord_elems(size_elem_list,dimension);
              IntList elems;
              // Fill a list elems and coordinates of the elems
              for (int ind_elem=0; ind_elem<size_elem_list; ind_elem++)
                {
                  int elem = elem_list_par_pos[ind_elem];
                  elems.add_if_not(elem);
                  for (int dir=0; dir<dimension; dir++)
                    coord_elems(ind_elem,dir) = xp(elem,dir);
                }
              int elem_identifie;
              if (size_elem_list==1)
                {
                  // Easy, one element in the list
                  elem_identifie = elems[0];
                }
              else if (size_elem_list==2 && norme_array(delt_dist)==0)
                {
                  // Check first if it is not 2 elements surrounding 2 internal superposed boundary faces, in this
                  // case, discards the element linked to the local boundary face
                  int local_face = fr_vf2.num_premiere_face()+ind_face;
                  int local_elem = zvf.face_voisins(local_face,0);
                  if (local_elem==-1) local_elem = zvf.face_voisins(local_face,1);
                  //Cerr << "[" << Process::me() << "<->" << pe << "] nb_faces=" << fr_vf2.nb_faces() << " face= " << ind_face << " " << local_elem << " " << elems[0] << " " << elems[1] << finl;
                  if (elems[0]==local_elem)
                    elem_identifie=elems[1]; // Remote elem
                  else if (elems[1]==local_elem)
                    elem_identifie=elems[0]; // Remote elem
                  else
                    {
                      // Surely 2 elements surrounding an internal face so select one element:
                      elem_identifie = Domaine::identifie_item_unique(elems,coord_elems,remote_point);
                    }
                }
              else
                {
                  // Find a unic remote elem:
                  elem_identifie = Domaine::identifie_item_unique(elems,coord_elems,remote_point);
                }
              elem_list[ind_face] = (elem_identifie < nb_elem_domaine1 ? elem_identifie : -1);
            }
        }
 
      // We fill these three arrays:
      DoubleTab& coord1 = inconnues1_coords_to_eval_[pe];
 
      ArrOfInt& num_elem = inconnues1_elems_[pe];
 
      ArrOfInt& index_to_recv = indexes_to_recv[pe];
 
 
      int nb_remote_faces = 0;
      // Loop on local faces on the process pe:
      for (int face = 0; face < nb_faces_on_pe; face++)
        {
          const int elem = elem_list[face];
          if (elem < 0 || elem >= nb_elem_domaine1)
            {
              // This coordinate is not on this processor...
            }
          else
            {
              coord1.resize(nb_remote_faces + 1, dim);
              for (int j = 0; j < dim; j++)
                coord1(nb_remote_faces, j) = remote_coords_on_pe(face, j);
              nb_remote_faces++;
              num_elem.append_array(elem);
              index_to_recv.append_array(face);
            }
        }
    }
 
  // Send index_to_recv data to each processor
  envoyer_all_to_all(indexes_to_recv, inconnues2_faces_);
 
  // Check that each coordinate is received from at least one processor:
  const int nb_faces2 = remote_coords[moi].dimension(0);
  ArrOfDouble count(nb_faces2);
  for (int pe = 0; pe < nprocs; pe++)
    {
      const ArrOfInt index_to_recv = inconnues2_faces_[pe];
      const int n = index_to_recv.size_array();
      for (int i = 0; i < n; i++)
        count[index_to_recv[i]] ++;
      //Cerr << index_to_recv << finl;
    }
  bool error_1 = false, error_2 = false;
  for (int i = 0; i < nb_faces2; i++)
    {
      if (count[i] < 1)
        {
          Journal() << "Champ_front_recyclage : error face (number "<< i<< ") " << remote_coords[moi](i, 0) << " " << remote_coords[moi](i, 1) << " "
                    << ((dim==3)?remote_coords[moi](i, 2):0.) << " no matching element found" << finl;
          error_1 = true;
        }
      else if (count[i]>1)
        {
          Journal() << "Champ_front_recyclage : error face (number "<< i<< ") " << remote_coords[moi](i, 0) << " " << remote_coords[moi](i, 1) << " "
                    << ((dim==3)?remote_coords[moi](i, 2):0.) << " " << count[i] << " corresponding points" << finl;
          error_2 = true;
        }
 
    }
  if (mp_or(error_1))
    {
      Cerr << "Error in Champ_front_recyclage.\n"
           << "No matching source element found for some faces of boundary " << nom_bord2
           << " (see log files)." << finl;
      Cerr << "If you are sure that boundary faces match with faces of the other domain, try adding this keyword:" << finl;
      Cerr << "PrecisionGeom 1.e-6" << finl;
      Cerr << "Just after the following keyword in your data file:" << finl;
      Cerr << "dimension " << Objet_U::dimension << finl;
      barrier();
      exit();
    }
 
  if (mp_or(error_2))
    {
      Cerr << "Error in Champ_front_recyclage.\n"
           << "At least two elements match for some faces of boundary " << nom_bord2
           << " (see log files)." << finl;
 
      barrier();
      exit();
    }
 
  initialiser_moyenne_imposee(moyenne_imposee_);
  initialiser_moyenne_recyclee(moyenne_recyclee_);
  mettre_a_jour(temps);
 
  return 1;
 
}
 
//-Evaluation of the values of l_inconnue1 at the points located
// in plane1 (xv(faces_bord2) + delt_dist)
//-Exchange of data between processors
//-Computation of moyenne_imposee_ and moyenne_recyclee_ if time-dependent evaluation
// (done by the local processor on bord2 from data received from other
//  processors or specified in the data set)
//
void Champ_front_recyclage::mettre_a_jour(double temps)
{
 
  DoubleTab& tab = valeurs_au_temps(temps);
  const int nprocs = nproc();
  DoubleTabs values_to_send(nprocs);
  DoubleTabs values_to_recv(nprocs);
  // temporary array (because valeur_aux_elems wants intvect and not arrofint)
  IntVect elems;
 
 
  for (int pe = 0; pe < nprocs; pe++)
    {
      elems.resize(inconnues1_elems_[pe].size_array());
      elems.inject_array(inconnues1_elems_[pe]);
 
      const int n = elems.size_array();
      values_to_send[pe].resize(n,nb_compo_);
      l_inconnue1->valeur_aux_elems(inconnues1_coords_to_eval_[pe],
                                    elems,
                                    values_to_send[pe]);
 
    }
  envoyer_all_to_all(values_to_send, values_to_recv);
  for (int pe = 0; pe < nprocs; pe++)
    {
      const DoubleTab& values = values_to_recv[pe];
      const ArrOfInt index_to_recv = inconnues2_faces_[pe];
      const int n = index_to_recv.size_array();
      for (int i = 0; i < n; i++)
        {
          const int idx = index_to_recv[i];
          for (int j = 0; j < nb_compo_; j++)
            tab(idx, j) = values(i, j);
        }
    }
 
  const Front_VF& fr_vf2 = ref_cast(Front_VF,la_frontiere_dis.valeur());
  int nb_faces2 = fr_vf2.nb_faces();
 
  calcul_moyenne_imposee(tab,temps);
  calcul_moyenne_recyclee(tab,temps);
 
  for (int i=0; i<nb_faces2; i++)
    for (int dir=0; dir<nb_compo_; dir++)
      tab(i,dir) = ampli_moy_imposee_(dir)*moyenne_imposee_(i,dir) + ampli_fluct_(dir)*(tab(i,dir)-ampli_moy_recyclee_(dir)*moyenne_recyclee_(i,dir));
 
  tab.echange_espace_virtuel();
 
}
 
//Building data necessary for estimating moyenne_imposee_
//according to the chosen evaluation method (methode_moy_impos_)
//
void Champ_front_recyclage::initialiser_moyenne_imposee(DoubleTab& moyenne)
{
 
  const Front_VF& fr_vf2 = ref_cast(Front_VF,la_frontiere_dis.valeur());
 
  if (methode_moy_impos_==1)
    {
      profil_2.dimensionner(nb_compo_);
      for (int i=0; i<nb_compo_; i++)
        {
          profil_2[i].setNbVar(4);
          profil_2[i].setString(fcts_xyz[i]);
          profil_2[i].addVar("t");
          profil_2[i].addVar("x");
          profil_2[i].addVar("y");
          profil_2[i].addVar("z");
          profil_2[i].parseString();
        }
    }
  else if (methode_moy_impos_==2)
    {
      lire_fichier_format1(moyenne_imposee_,fr_vf2,fich_impos_);
    }
  else if (methode_moy_impos_==3)
    {
      lire_fichier_format2(moyenne_imposee_,fr_vf2,fich_impos_);
    }
  else if (methode_moy_impos_==4)
    {
      lire_fichier_format3(moyenne_imposee_,fr_vf2,fich_impos_,fich_maillage_);
    }
  else if (methode_moy_impos_==5)
    {
      const Domaine_dis_base& domaine_dis2 = fr_vf2.domaine_dis();
      const Domaine_VF& le_dom2        = ref_cast(Domaine_VF,domaine_dis2);
      const DoubleTab&  xv2          = le_dom2.xv();
 
      double Kappa = 0.41;
      double d_paroi,y_plus,U_log;
 
      int nbface2 = fr_vf2.nb_faces();
      int num1 = fr_vf2.num_premiere_face();
      int num2 = num1 + nbface2;
      int num_face_loc;
 
      for (int num_face2=num1; num_face2<num2; num_face2++)
        {
          num_face_loc = num_face2-num1;
          d_paroi = diametre/2.;
 
          if(dimension==3)
            {
              if(abs(ndir)==1) d_paroi -= sqrt(xv2(num_face2,1)*xv2(num_face2,1) + xv2(num_face2,2)*xv2(num_face2,2));
              if(abs(ndir)==2) d_paroi -= sqrt(xv2(num_face2,0)*xv2(num_face2,0) + xv2(num_face2,2)*xv2(num_face2,2));
              if(abs(ndir)==3) d_paroi -= sqrt(xv2(num_face2,0)*xv2(num_face2,0) + xv2(num_face2,1)*xv2(num_face2,1));
            }
          else
            {
              if(abs(ndir)==1) d_paroi -= sqrt(xv2(num_face2,1)*xv2(num_face2,1));
              if(abs(ndir)==2) d_paroi -= sqrt(xv2(num_face2,0)*xv2(num_face2,0));
            }
 
          y_plus = d_paroi*u_tau/visco_cin;
          U_log  = (1./Kappa)*log(y_plus) + 5.1;
          U_log *=  u_tau;
          U_log *= ndir/abs(ndir);
 
          for (int dir=0; dir<nb_compo_; dir++)
            moyenne(num_face_loc,dir) = ((dir+1)==abs(ndir))*U_log;
        }
    }
}
 
//Building data necessary for estimating moyenne_recyclee_
//according to the chosen evaluation method (methode_moy_recycl_)
//
void Champ_front_recyclage::initialiser_moyenne_recyclee(DoubleTab& moyenne)
{
  const Front_VF& fr_vf2 = ref_cast(Front_VF,la_frontiere_dis.valeur());
  if (methode_moy_recycl_==2)
    {
      lire_fichier_format1(moyenne_recyclee_,fr_vf2,fich_recycl_);
    }
  else if (methode_moy_recycl_==3)
    {
      lire_fichier_format2(moyenne_recyclee_,fr_vf2,fich_recycl_);
    }
  else if (methode_moy_recycl_==4)
    {
      lire_fichier_format3(moyenne_recyclee_,fr_vf2,fich_recycl_,fich_maillage_);
    }
}
 
 
//Evaluation of moyenne_imposee_ according to the chosen evaluation method
//-methode_moy_impos_==1 : analytical profile
//-methode_moy_impos_==2 : no update (see initialiser_moyenne_imposee())
//                           read at initialization and interpolation of data read
//-methode_moy_impos_==3 : no update (see initialiser_moyenne_imposee())
//                           read at initialization to perform connexion_approchee
//-methode_moy_impos_==4 : no update (see initialiser_moyenne_imposee())
//                         read at initialization to perform connexion_exacte
//-methode_moy_impos_==5 : no update (see initialiser_moyenne_imposee())
//                           built at initialization using a log law:
//                           U_log  = ndir/abs(ndir)*u_tau*(1./Kappa)*log(y_plus) + 5.1;
//
void Champ_front_recyclage::calcul_moyenne_imposee(const DoubleTab& tab,double temps)
{
  if (methode_moy_impos_==1)
    {
      const Front_VF& fr_vf2 = ref_cast(Front_VF,la_frontiere_dis.valeur());
      const Domaine_dis_base& domaine_dis2 = fr_vf2.domaine_dis();
      const Domaine_VF&  le_dom2 = ref_cast(Domaine_VF,domaine_dis2);
      const DoubleTab&  xv2 = le_dom2.xv();
      int nb_faces_bord2 = fr_vf2.nb_faces();
      int ndeb = fr_vf2.num_premiere_face();
 
      for (int i=0; i<nb_faces_bord2; i++)
        for (int j=0; j<nb_compo_; j++)
          {
            profil_2[j].setVar(0,temps);
            for (int d=0; d<dimension; d++)
              profil_2[j].setVar(d+1,xv2(ndeb+i,d));
            moyenne_imposee_(i,j) =  profil_2[j].eval();
          }
    }
  else if (methode_moy_impos_==2)
    {
    }
  else if (methode_moy_impos_==3)
    {
    }
  else if (methode_moy_impos_==4)
    {
    }
  else if (methode_moy_impos_==5)
    {
    }
  else if (methode_moy_impos_!=-1)
    {
      Cerr<<"The method used to calculate moyenne_imposee_ is not definded."<<finl;
      exit();
    }
}
 
//Evaluation of moyenne_recyclee_ according to the chosen evaluation method
//-methode_moy_recycl_==1 : surface average on bord2 of the values estimated
//                          by l_inconnue1 in plan1
//-methode_moy_recycl_==2 : no update (see initialiser_moyenne_recyclee)
//                            read at initialization and interpolation of data read
//-methode_moy_recycl_==3 : no update (see initialiser_moyenne_recyclee())
//                            read at initialization to perform connexion_approchee
//-methode_moy_recycl_==4 : no update (see initialiser_moyenne_recyclee())
//                          read at initialization to perform connexion_exacte
//
void Champ_front_recyclage::calcul_moyenne_recyclee(const DoubleTab& tab,double temps)
{
  if (methode_moy_recycl_==1)
    {
      const Front_VF& fr_vf2 = ref_cast(Front_VF,la_frontiere_dis.valeur());
      int nb_faces_bord2 = fr_vf2.nb_faces();
      DoubleVect face_surfaces;
      const Faces& les_faces_bord = fr_vf2.frontiere().faces();
      les_faces_bord.calculer_surfaces(face_surfaces);
 
      DoubleVect moyenne(nb_compo_);
      double somme_si = 0.;
      moyenne = 0.;
      for (int i=0; i<nb_faces_bord2; i++)
        {
          for (int j=0; j<nb_compo_; j++)
            moyenne(j) += tab(i,j)*face_surfaces(i);
 
          somme_si += face_surfaces(i);
        }
      somme_si = mp_sum(somme_si);
 
      for (int j=0; j<nb_compo_; j++)
        moyenne(j) = mp_sum(moyenne(j));
      moyenne /= somme_si;
 
      for (int i=0; i<nb_faces_bord2; i++)
        for (int j=0; j<nb_compo_; j++)
          moyenne_recyclee_(i,j) = moyenne(j);
    }
  else if (methode_moy_recycl_==2)
    {
    }
  else if (methode_moy_recycl_==3)
    {
    }
  else if (methode_moy_recycl_==4)
    {
    }
  else if (methode_moy_recycl_!=-1)
    {
      Cerr<<"The method used to calculate moyenne_recyclee_ is not defined."<<finl;
      exit();
    }
}
 
//-Reading the file to dimension y_Pb2 and U_Pb2
//-Using the data read from the file for the three points
// with coordinates closest to y in order to interpolate
// the associated variable
//
double Champ_front_recyclage::UPb(double y,Nom nom_fich)
{
 
  int nb_pts;
  DoubleVect y_Pb2;
  DoubleVect U_Pb2;
  int l=0;
 
  if (Process::je_suis_maitre())
    {
      EFichier fich_Pb2(nom_fich);
      while (!fich_Pb2.eof())
        {
          y_Pb2.resize(l+1);
          U_Pb2.resize(l+1);
          fich_Pb2 >> y_Pb2(l) >> U_Pb2(l);
          l++;
        }
      fich_Pb2.close();
    }
 
  envoyer_broadcast(l,0);
  if (!(Process::je_suis_maitre()))
    {
      y_Pb2.resize(l);
      U_Pb2.resize(l);
    }
  envoyer_broadcast(y_Pb2,0);
  envoyer_broadcast(U_Pb2,0);
 
  nb_pts=l-1;
  y_Pb2.resize(nb_pts);
  U_Pb2.resize(nb_pts);
 
  l=0;
  int l_deb;
 
  while ((y>y_Pb2(l)) && (l<nb_pts))
    {
      l++;
    }
 
  if (l<2)
    {
      l_deb=0;
    }
  else
    {
      if (l>(nb_pts-2))
        {
          l_deb=nb_pts-3;
        }
      else
        {
          if ( (y-y_Pb2(l-1)) > (y_Pb2(l)-y) )
            {
              l_deb=l-2;
            }
          else
            {
              l_deb=l-1;
            }
        }
    }
 
  double ya=y_Pb2(l_deb);
  double Ua=U_Pb2(l_deb);
  double yb=y_Pb2(l_deb+1);
  double Ub=U_Pb2(l_deb+1);
  double yc=y_Pb2(l_deb+2);
  double Uc=U_Pb2(l_deb+2);
  double U;
 
  U=(y-ya)*(y-yb)/(yc-ya)/(yc-yb)*Uc;
  U+=(y-yb)*(y-yc)/(ya-yb)/(ya-yc)*Ua;
  U+=(y-yc)*(y-ya)/(yb-yc)/(yb-ya)*Ub;
 
  return U;
 
}
 
//-Determine a direction orthogonal to the anisotropic direction
//-For each evaluation point located in plane1:
//  Estimate its coordinate y in the orthogonal direction
//  Initialize the mean (in the anisotropic direction) from
//  the UPb function for coordinate y
//
void Champ_front_recyclage::lire_fichier_format1(DoubleTab& moyenne,
                                                 const Frontiere_dis_base& fr_vf,
                                                 const Nom& nom_fich)
{
  const Front_VF& fr_vf2 = ref_cast(Front_VF,fr_vf);
  DoubleTab coords;
  get_coord_faces(fr_vf2,coords,delt_dist);
  const int nb_faces2 = coords.dimension(0);
 
  int dir_ortho = -1;
  if ((ndir==1) || (ndir==2))
    dir_ortho = ndir;
  else if (ndir==3)
    dir_ortho = 0;
  else
    {
      Cerr<<"For this method the anisotropic direction mut be 1 2 or 3"<<finl;
      exit();
    }
 
  for (int i=0; i<nb_faces2; i++)
    {
      double y = coords(i,dir_ortho);
      moyenne(i,ndir-1) = UPb(y,nom_fich);
    }
}
 
//-For each face of bord2 (coordinates x2 y2 z2):
//  Read from file nom_fich the coordinates of points (x1 y1 z1)
//  and the value(s) of the associated variable var1 [var2] [var3]
//  Initialize the mean (in the anisotropic direction if specified,
//  or for all components otherwise) from var1 [var2] [var3]
//  associated to point x1 y1 z1 which is at the smallest distance from x2 y2 z2
//
// Possible parallel issue if (x1 y1 z1) do not cover the virtual parts
// or if two points are at equal distance from a point (x2 y2 z2)
//
void Champ_front_recyclage::lire_fichier_format2(DoubleTab& moyenne,
                                                 const Frontiere_dis_base& fr_vf,
                                                 const Nom& nom_fich)
{
 
  int N1;
  int num_face1, indice_compo;
  int num_face2, num_face2_loc;
  double x1,y1,z1;
  double x2,y2,z2;
  z2 = 0.;
  DoubleVect val(nb_compo_);
 
  // Read the file and keep in memory
  ifstream fic_Utargtest(nom_fich); // test
  if (!fic_Utargtest)
    {
      Cerr <<"No file of name : "<<nom_fich<<" is available."<<finl;
      exit();
    }
  LecFicDiffuse fic_Utarg(nom_fich, ios::in);
  fic_Utarg >> N1;
  ArrOfDouble fich_buffer(N1*(nb_compo_+dimension));
  fic_Utarg.get(fich_buffer.addr(), fich_buffer.size_array());
  fic_Utarg.close();
 
  const Front_VF& fr_vf2 = ref_cast(Front_VF,fr_vf);
  int num1 = fr_vf2.num_premiere_face();
  int num2 = num1 + fr_vf2.nb_faces();
  const Domaine_dis_base& domaine_dis2 = fr_vf2.domaine_dis();
  const Domaine_VF& zvf2 = ref_cast(Domaine_VF,domaine_dis2);
  const DoubleTab& xv2 = zvf2.xv();
 
  for (num_face2=num1; num_face2<num2; num_face2++)
    {
      x2 = xv2(num_face2,0);
      y2 = xv2(num_face2,1);
      if (dimension==3) z2 = xv2(num_face2,2);
      num_face2_loc = num_face2-num1;
 
      double dist = 1.E6;
      double dist_eval;
 
      for (num_face1=0; num_face1<N1; num_face1++)
        {
          x1 = fich_buffer[num_face1*(nb_compo_+dimension)];
          y1 = fich_buffer[num_face1*(nb_compo_+dimension)+1];
          if(dimension==3) z1 = fich_buffer[num_face1*(nb_compo_+dimension)+2];
          dist_eval = (x2-x1)*(x2-x1)+(y2-y1)*(y2-y1);
          if (dimension==3) dist_eval+=(z2-z1)*(z2-z1);
 
          if(dist_eval<dist)
            {
              for (indice_compo=0; indice_compo<nb_compo_; indice_compo++)
                val(indice_compo) = fich_buffer[num_face1*(nb_compo_+dimension)+dimension+indice_compo];
              dist = dist_eval;
              if ((abs(ndir)<4) && (abs(ndir)>0))
                moyenne(num_face2_loc,abs(ndir)-1) = val(abs(ndir)-1);
              else
                for (indice_compo=0; indice_compo<nb_compo_; indice_compo++)
                  moyenne(num_face2_loc,indice_compo) = val(indice_compo);
            }
        }
    }
}
 
//-For each face of bord2 (coordinates x2 y2 z2):
//  Read from file nom_fich2 the coordinates of points (x1 y1 z1)
//  Read from file nom_fich1 the values of the associated variable
//  and initialize the mean from these values. Reading stops
//  when the values associated to point x1 y1 z1 that is in exact
//  connection with x2 y2 z2 along the anisotropic direction are read
//
void Champ_front_recyclage::lire_fichier_format3(DoubleTab& moyenne,
                                                 const Frontiere_dis_base& fr_vf,
                                                 const Nom& nom_fich1,const Nom& nom_fich2)
{
  int Nbfaces,trouve,compteur;
  int num_face2,num_face2_loc;
  double x1=0.,y1=0.,z1=0.;
  double x2=0.,y2=0.,z2=0.;
  double eps=1.e-4;
 
  Nom fichier = "VERIF_face_perio.";
  fichier += nom_pb1;
  EcrFicCollecte fic(fichier);
 
  ifstream fic_geom_test(nom_fich2);
  if (!fic_geom_test)
    {
      Cerr <<"No file of name : "<<nom_fich2<<" is available."<<finl;
      exit();
    }
 
  ifstream fic_Umoy_test(nom_fich1);
  if (!fic_Umoy_test)
    {
      Cerr <<"No file of name : "<<nom_fich1<<" is available."<<finl;
      exit();
    }
 
  if (abs(ndir)!=1 && abs(ndir)!=2 && abs(ndir)!=3)
    {
      Cerr << finl;
      Cerr << "Error while reading the boundary condition " << fr_vf.le_nom() << finl;
      Cerr << "You must define the keyword \"direction_anisotrope\"" << finl;
      exit();
    }
 
  if (abs(ndir)==3 && dimension==2)
    {
      Cerr << finl;
      Cerr << "Error while reading the boundary condition " << fr_vf.le_nom() << finl;
      Cerr << "You must define the keyword \"direction_anisotrope\" with a direction different than 'z' in 2D" << finl;
      exit();
    }
 
  const Front_VF& fr_vf2 = ref_cast(Front_VF,fr_vf);
  int num1 = fr_vf2.num_premiere_face();
  int num2 = num1 + fr_vf2.nb_faces();
  const Domaine_dis_base& domaine_dis2 = fr_vf2.domaine_dis();
  const Domaine_VF& zvf2 = ref_cast(Domaine_VF,domaine_dis2);
  const DoubleTab& xv2 = zvf2.xv();
 
  // Reading and buffering the two input files
  LecFicDiffuse fic_geom(nom_fich2);
  LecFicDiffuse fic_Umoy(nom_fich1);
  fic_geom >> Nbfaces;
  fic_Umoy >> Nbfaces; // useless but we must read it
  ArrOfDouble fich_geom_buffer(Nbfaces*(dimension+1));
  ArrOfDouble fich_Umoy_buffer(Nbfaces*(dimension+1));
  fic_geom.get(fich_geom_buffer.addr(), fich_geom_buffer.size_array());
  fic_Umoy.get(fich_Umoy_buffer.addr(), fich_Umoy_buffer.size_array());
  fic_geom.close();
  fic_Umoy.close();
 
  for (num_face2=num1; num_face2<num2; num_face2++)
    {
      x2 = xv2(num_face2,0);
      y2 = xv2(num_face2,1);
      if (dimension==3) z2 = xv2(num_face2,2);
      num_face2_loc = num_face2-num1;
 
      trouve = 0;
      compteur = 0;
 
      // Looking for (x2,y2[,z2]) into (x1,y1[,z1])
      while (compteur<Nbfaces)
        {
          x1 = fich_geom_buffer[compteur*(dimension+1)+1];
          y1 = fich_geom_buffer[compteur*(dimension+1)+2];
          if (dimension==3)
            {
              z1 = fich_geom_buffer[compteur*(dimension+1)+3];
              if(abs(ndir)==1 && (est_egal(y1,y2,eps)) && (est_egal(z1,z2,eps))) trouve=1;
              if(abs(ndir)==2 && (est_egal(x1,x2,eps)) && (est_egal(z1,z2,eps))) trouve=1;
              if(abs(ndir)==3 && (est_egal(x1,x2,eps)) && (est_egal(y1,y2,eps))) trouve=1;
            }
          else
            {
              if(abs(ndir)==1 && (est_egal(y1,y2,eps))) trouve=1;
              if(abs(ndir)==2 && (est_egal(x1,x2,eps))) trouve=1;
            }
          if (trouve)
            {
              moyenne(num_face2_loc,0) = fich_Umoy_buffer[compteur*(dimension+1)+1];
              moyenne(num_face2_loc,1) = fich_Umoy_buffer[compteur*(dimension+1)+2];
              if (dimension==3) moyenne(num_face2_loc,2) = fich_Umoy_buffer[compteur*(dimension+1)+3];
              break;
            }
          compteur++;
        }
 
      if (trouve==0)
        {
          Cerr << " the face "<<num_face2<<" with coordinates : x="<<x2<<" y="<<y2;
          if (dimension==3)
            Cerr << " z="<<z2;
          Cerr << "\n has no corresponding faces in "<<nom_fich2<<finl;
          exit();
        }
 
      if (dimension==3)
        {
          fic <<" num_face1 " << fich_geom_buffer[compteur*(dimension+1)]  << " " << x1 << " " << y1 << " " << z1 << finl;
          fic <<" num_face2 " << num_face2  << " " << x2 << " " << y2 << " " << z2 << finl;
          fic <<" "  << finl;
        }
      else
        {
          fic <<" num_face1 " << fich_geom_buffer[compteur*(dimension+1)]  << " " << x1 << " " << y1  << finl;
          fic <<" num_face2 " << num_face2  << " " << x2 << " " << y2  << finl;
          fic <<" "  << finl;
        }
 
    }
  fic.close();
}