Champ_Face_PolyMAC_MPFA.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_Face_PolyMAC_MPFA.h>
#include <Domaine.h>
#include <Domaine_VF.h>
#include <Champ_Uniforme.h>
#include <Domaine_PolyMAC_MPFA.h>
#include <Domaine_Cl_PolyMAC_family.h>
#include <TRUSTLists.h>
#include <Dirichlet.h>
#include <Dirichlet_homogene.h>
#include <Symetrie.h>
#include <EChaine.h>
#include <TRUSTTab_parts.h>
#include <Linear_algebra_tools_impl.h>
#include <Probleme_base.h>
#include <Equation_base.h>
#include <Schema_Temps_base.h>
#include <Champ_Fonc_reprise.h>
#include <array>
#include <cmath>
#include <Pb_Multiphase.h>
 
namespace
{
inline void project_axis_face(const IntTab& f_s, const DoubleTab& xs, const DoubleTab& xp,
                              DoubleTrav& xfb, int f, int e)
{
  int s0 = -1, s1 = -1;
  for (int i = 0; i < f_s.dimension(1) && s1 < 0; i++)
    {
      const int s = f_s(f, i);
      if (s < 0) break;
      if (s0 < 0) s0 = s;
      else s1 = s;
    }
 
  assert(s0 >= 0 && s1 >= 0);
  const double tx = xs(s1, 0) - xs(s0, 0), ty = xs(s1, 1) - xs(s0, 1);
  const double nx = ty, ny = -tx;
  const double n2 = nx * nx + ny * ny;
 
  const double mx = 0.5 * (xs(s0, 0) + xs(s1, 0)), my = 0.5 * (xs(s0, 1) + xs(s1, 1));
  const double dx = xp(e, 0) - mx, dy = xp(e, 1) - my;
  const double coeff = (dx * nx + dy * ny) / n2;
 
  xfb(f, 0) = xp(e, 0) - coeff * nx;
  xfb(f, 1) = xp(e, 1) - coeff * ny;
}
}
 
Implemente_instanciable(Champ_Face_PolyMAC_MPFA,"Champ_Face_PolyMAC_MPFA",Champ_Face_PolyMAC_HFV) ;
 
Sortie& Champ_Face_PolyMAC_MPFA::printOn(Sortie& os) const { return os << que_suis_je() << " " << le_nom(); }
 
Entree& Champ_Face_PolyMAC_MPFA::readOn(Entree& is) { return is; }
 
int Champ_Face_PolyMAC_MPFA::fixer_nb_valeurs_nodales(int n)
{
  // use the same kind of code as in Champ_Fonc_P0_base, but retrieve the number of faces instead of the number of elements
  // note: the code used in Champ_Fonc_P0_base::fixer_nb_valeurs_nodales() to retrieve the discrete domain is somewhat surprising...
 
  assert(n == domaine_PolyMAC_MPFA().nb_faces() || n < 0);
 
  // Problem: nb_comp equals dimension but we only want one dimension !!!
  // HACK:
  int old_nb_compo = nb_compo_;
  if(nb_compo_ != 1) nb_compo_ /= dimension;
 
  /* variables: normal values at faces, then values at elements in blocks -> so that line_size() works */
  creer_tableau_distribue(domaine_PolyMAC_MPFA().md_vector_faces());
  nb_compo_ = old_nb_compo;
  return n;
}
 
void Champ_Face_PolyMAC_MPFA::init_auxiliary_variables()
{
  for (int n = 0; n < nb_valeurs_temporelles(); n++)
    {
      DoubleTab& vals = futur(n);
      vals.set_md_vector(MD_Vector()); //remove the MD_Vector...
      vals.resize_dim0(domaine_PolyMAC_MPFA().mdv_ch_face->get_nb_items_tot()); //...resize to the correct size...
      vals.set_md_vector(domaine_PolyMAC_MPFA().mdv_ch_face); //...and restore the correct MD_Vector
      update_ve(vals);
    }
}
 
int Champ_Face_PolyMAC_MPFA::reprendre(Entree& fich)
{
  if (! via_ch_fonc_reprise()) return Champ_Inc_base::reprendre(fich); /* ie: resume last time ! */
 
  const Pb_Multiphase * pbm = mon_equation_non_nul() ? (sub_type(Pb_Multiphase, equation().probleme()) ? &ref_cast(Pb_Multiphase, equation().probleme()) : nullptr) : nullptr;
  if (pbm) return Champ_Inc_base::reprendre(fich);
 
  // otherwise do this ...
  const Domaine_PolyMAC_MPFA* domaine = le_dom_VF ? &ref_cast( Domaine_PolyMAC_MPFA,le_dom_VF.valeur()) : nullptr;
  valeurs().set_md_vector(MD_Vector()); //remove the MD_Vector...
  valeurs().resize(0);
  int ret = Champ_Inc_base::reprendre(fich);
  //restore the correct one if possible
  if (domaine) valeurs().set_md_vector(valeurs().dimension_tot(0) > domaine->nb_faces_tot() ? domaine->mdv_faces_aretes : domaine->md_vector_faces());
  return ret;
}
 
Champ_base& Champ_Face_PolyMAC_MPFA::affecter_(const Champ_base& ch)
{
  const DoubleTab& ch_val = ch.valeurs();
  DoubleTab& val = valeurs();
  const int N = val.line_size(), D = Objet_U::dimension;
 
  if (sub_type(Champ_Fonc_reprise, ch))
    {
      if (val.dimension_tot(0) == ch_val.dimension_tot(0) && val.line_size() == ch_val.line_size())
        val = ch_val;
      else
        {
          init_auxiliary_variables();
          for (int i = 0; i < ch_val.dimension_tot(0); i++)
            for (int j = 0; j < ch_val.line_size(); j++)
              val(i,j) = ch_val(i,j);
        }
    }
  else
    {
      const Domaine_PolyMAC_MPFA& domaine = domaine_PolyMAC_MPFA();
      const DoubleVect& fs = domaine.face_surfaces();
      const DoubleTab& nf = domaine.face_normales(), &xv = domaine.xv();
      const int unif = sub_type(Champ_Uniforme, ch);
      DoubleTab eval;
 
      if (unif)
        eval = ch_val;
      else
        eval.resize(val.dimension_tot(0), N * D), ch.valeur_aux(xv, eval);
 
      for (int f = 0; f < domaine.nb_faces_tot(); f++)
        for (int d = 0; d < D; d++)
          for (int n = 0; n < N; n++)
            if (fs(f) > 0)
              val(f, n) += eval(unif ? 0 : f, N * d + n) * nf(f, d) / fs(f);
 
      update_ve(val);
      //copy into all time slots
      val.echange_espace_virtuel();
      for (int i = 1; i < les_valeurs->nb_cases(); i++)
        les_valeurs[i].valeurs() = val;
    }
  return *this;
}
 
/**
 * @brief Interpolates face-based velocity to element centers using a first-order method.
 *
 * Computes the element-based velocity components by performing a weighted average
 * of neighboring face velocities, projected onto the vector from face center to element center.
 *
 * This method uses a first-order accurate stencil and does not include boundary corrections.
 *
 * @param[out] val Output field containing element-centered velocity components
 *                 in the last D * number_of_elements entries.
 */
void Champ_Face_PolyMAC_MPFA::update_ve(DoubleTab& val) const
{
  const Domaine_PolyMAC_MPFA& domaine = domaine_PolyMAC_MPFA();
  if (valeurs().get_md_vector() != domaine.mdv_ch_face)
    return; //no auxiliary variables -> nothing to do
 
  const DoubleVect& fs = domaine.face_surfaces(), &ve = domaine.volumes(),
                    *pf = mon_equation_non_nul() ? &equation().milieu().porosite_face() : nullptr,
                     *pe = pf ? &equation().milieu().porosite_elem() : nullptr;
 
  const DoubleTab& xp = domaine.xp(), &xv = domaine.xv();
  const IntTab& e_f = domaine.elem_faces(), &f_e = domaine.face_voisins();
  const int D = dimension, N = val.line_size(), ne_tot = domaine.nb_elem_tot(), nf_tot = domaine.nb_faces_tot();
 
  double fac;
  for (int e = 0; e < ne_tot; e++)
    {
      for (int d = 0; d < D; d++)
        for (int n = 0; n < N; n++)
          val(nf_tot + D * e + d, n) = 0;
      for (int j = 0; j < e_f.dimension(1); j++)
        {
          const int f = e_f(e, j);
          if (f < 0) continue;
 
          fac = (e == f_e(f, 0) ? 1 : -1) * (pf ? (*pf)(f) : 1.0) * fs(f) / ((pe ? (*pe)(e) : 1.0) * ve(e));
          for (int d = 0; d < D; d++)
            for (int n = 0; n < N; n++)
              val(nf_tot + D * e + d, n) += fac * (xv(f, d) - xp(e, d)) * val(f, n);
        }
    }
}
 
/**
 * @brief Initializes second-order interpolation coefficients from faces to element centers.
 *
 * This routine computes the interpolation stencil and weights for reconstructing
 * the velocity at element centers with second-order accuracy.
 *
 * - Builds a least-squares system to correct the first-order estimate.
 * - Handles internal and boundary faces (Neumann, Dirichlet).
 * - Uses extended connectivity (neighboring elements via shared vertices).
 *
 * Resulting coefficients are stored internally for later use in update_ve2().
 */
void Champ_Face_PolyMAC_MPFA::init_ve2() const
{
  const Domaine_PolyMAC_MPFA& domaine = domaine_PolyMAC_MPFA();
 
  if (ve2d.dimension(0) || valeurs().get_md_vector() != domaine.mdv_ch_face)
    return; //already initialized or no auxiliary variables
 
  const DoubleVect& pf = equation().milieu().porosite_face(),
                    &pe = equation().milieu().porosite_elem(),
                     &fs = domaine.face_surfaces(), &ve = domaine.volumes();
  const DoubleTab& xp = domaine.xp(), &xv = domaine.xv(), &nf = domaine.face_normales(),
                   &xs = domaine.domaine().coord_sommets();
  const IntTab& f_e = domaine.face_voisins(), &e_f = domaine.elem_faces(),
                &e_s = domaine.domaine().les_elems(), &f_s = domaine.face_sommets();
 
  int D = dimension, nl = D * (D + 1), infoo = 0, un = 1;
  double eps = 1e-8, fac;
  const double *xf;
 
  init_fcl();
 
  //position of boundary face points: centroid if internal or Dirichlet, projection if Neumann
  DoubleTrav xfb(domaine.nb_faces_tot(), D), ve2, ve2i, A, B, P, W(1);
  IntTrav pvt;
 
  const double tol_norm2 = 1e-24;
  for (int f = 0; f < domaine.nb_faces_tot(); f++)
    if (fcl_(f, 0) == 1 || fcl_(f, 0) == 2) //Neumann / Symetrie
      {
        const int e = f_e(f, 0);
        const double nf_norm2 = domaine.dot(&nf(f, 0), &nf(f, 0));
 
        if (!bidim_axi || nf_norm2 > tol_norm2)
          {
            double scal = domaine.dot(&xv(f, 0), &nf(f, 0), &xp(e, 0)) / (fs(f) * fs(f));
            for (int d = 0; d < D; d++)
              xfb(f, d) = xp(e, d) + scal * nf(f, d);
          }
        else
          project_axis_face(f_s, xs, xp, xfb, f, e);
      }
    else if (fcl_(f, 0))
      {
        for (int d = 0; d < D; d++)
          xfb(f, d) = xv(f, d); //Dirichlet
      }
 
  /* vertex-element and element-element connectivities */
  std::vector<std::set<int>> s_f(domaine.domaine().nb_som()), e_s_f(domaine.nb_elem());
  for (int f = 0; f < domaine.nb_faces_tot(); f++)
    for (int i = 0; i < f_s.dimension(1); i++)
      {
        const int s = f_s(f, i);
        if (s < 0) continue;
 
        if (s < domaine.domaine().nb_som())
          s_f[s].insert(f);
      }
 
  for (int e = 0; e < domaine.nb_elem(); e++)
    for (int i = 0; i < e_s.dimension(1); i++)
      {
        const int s = e_s(e, i);
        if (s < 0) continue;
 
        for (auto &&fa : s_f[s])
          e_s_f[e].insert(fa);
      }
 
  ve2d.resize(1, 2);
  ve2bj.resize(0, 2);
  std::map<std::array<int, 2>, int> v_i; // v_i[{f, -1 (internal) or component }] = index
  std::vector<std::array<int, 2>> i_v; // v_i[i_v[f]] = f
 
  for (int e = 0; e < domaine.nb_elem(); e++, v_i.clear(), i_v.clear())
    {
      /* stencil: faces of the element and its vertex-connected neighbors + all components at its boundary faces */
      for (auto &&fa : e_s_f[e])
        if (!v_i.count( { { fa, -1 } }))
      v_i[ { { fa, -1 } }] = (int) i_v.size(), i_v.push_back( { { fa, -1 } });
      for (int i = 0; i < e_f.dimension(1) ; i++)
        {
          const int f = e_f(e, i);
          if (f < 0) continue;
 
          if (fcl_(f, 0))
            for (int d = 0; d < D; d++)
              v_i[ { { f, d } }] = (int) i_v.size(), i_v.push_back( { { f, d } });
        }
 
      /* first-order interpolation coefficients and weights (as in Domaine_PolyMAC_MPFA::{e,f}grad) */
      const int nc = (int) i_v.size();
 
      ve2.resize(nc, D);
      A.resize(nc, nl);
      B.resize(nc);
      P.resize(nc);
      pvt.resize(nc);
      ve2 = 0.;
 
      for (int i = 0; i < e_f.dimension(1); i++)
        {
          const int f = e_f(e, i);
          if (f < 0) continue;
 
          fac = (e == f_e(f, 0) ? 1 : -1) * fs(f) / ve(e);
          int j = v_i.at( { { f, -1 } });
 
          for (int d = 0; d < D; d++)
            ve2(j, d) += fac * (xv(f, d) - xp(e, d));
        }
 
      for (int i = 0; i < nc; i++)
        {
          int f = i_v[i][0];
          xf = i_v[i][1] < 0 ? &xv(f, 0) : &xfb(f, 0);
          P(i) = 1. / sqrt(domaine.dot(xf, xf, &xp(e, 0), &xp(e, 0)));
        }
 
      /* per component: correction to achieve second-order accuracy */
      for (int d = 0; d < D; d++)
        {
          /* systeme A.x = b */
          B = 0;
          pvt = 0;
 
          for (int i = 0; i < nc; i++)
            {
              int f = i_v[i][0];
              if (std::fabs(fs(f)) < 1e-20) continue;
              int db = i_v[i][1];
              xf = db < 0 ? &xv(f, 0) : &xfb(f, 0);
 
              int jb = 0;
              for (int j = 0; j < D; j++)
                for (int k = 0; k <= D; k++, jb++)
                  {
                    fac = (db < 0 ? nf(f, j) / fs(f) : (db == j)) * (k < D ? xf[k] - xp(e, k) : 1);
                    A(i, jb) = fac * P(i);
                    if (k < D)
                      B(jb) -= fac * ve2(i, d); //first-order interpolation error to be corrected
                  }
            }
 
          /* x de norme L2 minimale par dgels */
          int nw = -1, rank=-1;
 
          F77NAME(dgelsy)(&nl, &nc, &un, &A(0, 0), &nl, &B(0), &nc, &pvt(0), &eps, &rank, &W(0), &nw, &infoo);
 
          W.resize(nw = (int) std::lrint(W(0)));
 
          F77NAME(dgelsy)(&nl, &nc, &un, &A(0, 0), &nl, &B(0), &nc, &pvt(0), &eps, &rank, &W(0), &nw, &infoo);
          assert(infoo == 0);
 
          /* ajout dans ve2 */
          for (int i = 0; i < nc; i++)
            ve2(i, d) += P(i) * B(i);
        }
 
      /* implicit treatment of Neumann / Symmetry boundary conditions */
      Matrice33 M(1, 0, 0, 0, 1, 0, 0, 0, 1), iM;
 
      for (int i = 0; i < nc; i++)
        if (i_v[i][1] >= 0 && fcl_(i_v[i][0], 0) < 2)
          for (int j = 0; j < D; j++)
            M(j, i_v[i][1]) -= ve2(i, j);
 
      Matrice33::inverse(M, iM);
      ve2i.resize(nc, D);
 
      for (int i = 0; i < nc; i++)
        for (int j = 0; j < D; j++)
          {
            ve2i(i, j) = 0;
            for (int k = 0; k < D; k++)
              ve2i(i, j) += iM(j, k) * ve2(i, k);
          }
 
      /* storage */
      for (int d = 0; d < D; d++, ve2d.append_line(ve2c.size(), ve2bc.size()))
        for (int i = 0; i < nc; i++)
          if (std::fabs(ve2i(i, d)) > 1e-6 && (i_v[i][1] < 0 || fcl_(i_v[i][0], 0) == 3))
            {
              i_v[i][1] < 0 ? ve2j.append_line(i_v[i][0]) : ve2bj.append_line(i_v[i][0], i_v[i][1]);
              (i_v[i][1] < 0 ? &ve2c : &ve2bc)->append_line(ve2i(i, d) * pf(i_v[i][0]) / pe(e));
            }
    }
 
  CRIMP(ve2d);
  CRIMP(ve2j);
  CRIMP(ve2bj);
  CRIMP(ve2c);
  CRIMP(ve2bc);
}
 
/**
 * @brief Applies the second-order interpolation from face velocities to element centers.
 *
 * Uses coefficients computed by init_ve2() to interpolate the face-based velocity field
 * into element-centered values. Supports both internal faces and Dirichlet boundary conditions.
 *
 * @param[out] val Interpolated velocity at element centers.
 *                 The last D * number_of_elements entries are updated.
 * @param[in] incr If non-zero, Dirichlet boundary contributions are ignored (for incremental updates).
 */
void Champ_Face_PolyMAC_MPFA::update_ve2(DoubleTab& val, int incr) const
{
  const Domaine_PolyMAC_MPFA& domaine = domaine_PolyMAC_MPFA();
  if (valeurs().get_md_vector() != domaine.mdv_ch_face)
    return; //no auxiliary variables -> exit
 
  const Conds_lim& cls = domaine_Cl_dis().les_conditions_limites();
  const int D = dimension,  N = val.line_size(), nf_tot = domaine.nb_faces_tot();
 
  init_ve2();
  init_fcl();
 
  int ed = 0, i = nf_tot;
 
  for (int e = 0 ; e < domaine.nb_elem(); e++)
    for (int d = 0; d < D; d++, ed++, i++)
      for (int n = 0; n < N; n++)
        {
          /* "internal" part */
          val(i, n) = 0;
 
          for (int j = ve2d(ed, 0); j < ve2d(ed + 1, 0); j++)
            val(i, n) += ve2c(j) * val(ve2j(j), n);
 
          /* "Dirichlet boundary faces" part (except when computing increments) */
          if (!incr)
            for (int j = ve2d(ed, 1); j < ve2d(ed + 1, 1); j++)
              if (sub_type(Dirichlet, cls[fcl_(ve2bj(j, 0), 1)].valeur()))
                val(i, n) += ve2bc(j) * ref_cast(Dirichlet, cls[fcl_(ve2bj(j, 0), 1)].valeur()).val_imp(fcl_(ve2bj(j, 0), 2), N * ve2bj(j, 1) + n);
        }
 
  val.echange_espace_virtuel();
}
 
DoubleTab& Champ_Face_PolyMAC_MPFA::valeur_aux_elems(const DoubleTab& positions, const IntVect& les_polys, DoubleTab& val_elem) const
{
  if (valeurs().get_md_vector() != domaine_PolyMAC_MPFA().mdv_ch_face)
    return Champ_Face_PolyMAC_HFV::valeur_aux_elems_(valeurs(), positions, les_polys, val_elem);
  else
    return valeur_aux_elems_(le_champ().valeurs(), positions, les_polys, val_elem);
}
 
DoubleTab& Champ_Face_PolyMAC_MPFA::valeur_aux_elems_passe(const DoubleTab& positions, const IntVect& les_polys, DoubleTab& val_elem) const
{
  if (valeurs().get_md_vector() != domaine_PolyMAC_MPFA().mdv_ch_face)
    return Champ_Face_PolyMAC_HFV::valeur_aux_elems_(passe(), positions, les_polys, val_elem);
  else
    return valeur_aux_elems_(le_champ().passe(), positions, les_polys, val_elem);
}
 
DoubleTab& Champ_Face_PolyMAC_MPFA::valeur_aux_elems_(const DoubleTab& val_face, const DoubleTab& positions, const IntVect& les_polys, DoubleTab& val_elem) const
{
  const Domaine_PolyMAC_MPFA& domaine = domaine_PolyMAC_MPFA();
  int nb_compo = le_champ().nb_comp(), nf_tot = domaine.nb_faces_tot(), d, D = dimension, n, N = val_face.line_size();
 
  assert((positions.dimension(0) == les_polys.size()) || (positions.dimension_tot(0) == les_polys.size()));
 
  if (val_elem.nb_dim() > 2)
    Process::exit("TRUST error in Champ_Face_PolyMAC_MPFA::valeur_aux_elems_ : The DoubleTab val has more than 2 entries !");
 
  if (nb_compo == 1)
    Process::exit("TRUST error in Champ_Face_PolyMAC_MPFA::valeur_aux_elems_ : A scalar field cannot be of Champ_Face type !");
 
  for (int p = 0, e; p < les_polys.size(); p++)
    for (e = les_polys(p), d = 0; e < domaine.nb_elem() && d < D; d++)
      for (n = 0; n < N; n++)
        val_elem(p, N * d + n) = val_face(nf_tot + D * e + d, n);
  return val_elem;
}
 
DoubleVect& Champ_Face_PolyMAC_MPFA::valeur_aux_elems_compo(const DoubleTab& positions, const IntVect& polys, DoubleVect& val, int ncomp) const
{
  if (valeurs().get_md_vector() != domaine_PolyMAC_MPFA().mdv_ch_face)
    return Champ_Face_PolyMAC_HFV::valeur_aux_elems_compo(positions, polys, val, ncomp);
  int nf_tot = domaine_PolyMAC_MPFA().nb_faces_tot(), D = dimension, N = valeurs().line_size();
  assert(val.size_totale() >= polys.size());
 
  DoubleTab vfe(valeurs());
  update_ve(vfe);
 
  for (int p = 0, e; p < polys.size(); p++)
    val(p) = (e = polys(p)) < 0 ? 0. : vfe.addr()[N * (nf_tot + D * e) + ncomp];
 
  return val;
}
 
DoubleTab& Champ_Face_PolyMAC_MPFA::trace(const Frontiere_dis_base& fr, DoubleTab& x, double t, int distant) const
{
  assert(distant == 0);
  const Domaine_PolyMAC_MPFA& domaine = domaine_PolyMAC_MPFA();
  if (valeurs().get_md_vector() != domaine.mdv_ch_face)
    return Champ_Face_PolyMAC_HFV::trace(fr, x, t, distant);
 
  const bool vectoriel = (le_champ().nb_comp() > 1);
  const int dim = vectoriel ? dimension : 1, ne_tot = domaine.nb_elem_tot(), nf_tot = domaine.nb_faces_tot();
  const Front_VF& fr_vf = ref_cast(Front_VF, fr);
  const IntTab& face_voisins = domaine.face_voisins();
  const DoubleTab& val = valeurs(t);
 
  for (int i = 0; i < fr_vf.nb_faces(); i++)
    {
      const int face = fr_vf.num_premiere_face() + i, elem = face_voisins(face, 0);
      for (int d = 0; d < dim; d++)
        x(i, d) = val(vectoriel ? nf_tot + ne_tot * d + elem : face);
    }
 
  return x;
}