Op_Conv_EF_Stab_PolyMAC_HFV_Elem.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 <Op_Conv_EF_Stab_PolyMAC_HFV_Elem.h>
#include <Champ_Elem_PolyMAC_HFV.h>
#include <Champ_Face_PolyMAC_HFV.h>
#include <Masse_PolyMAC_HFV_Elem.h>
#include <Convection_Diffusion_std.h>
#include <Discretisation_base.h>
#include <Dirichlet_homogene.h>
#include <Domaine_Cl_PolyMAC_family.h>
#include <Schema_Temps_base.h>
#include <Domaine_Poly_base.h>
#include <Option_PolyMAC_family.h>
#include <Pb_Multiphase.h>
#include <Probleme_base.h>
#include <Matrix_tools.h>
#include <Milieu_base.h>
#include <Array_tools.h>
#include <TRUSTLists.h>
#include <Dirichlet.h>
#include <Symetrie.h>
#include <Param.h>
#include <cfloat>
#include <vector>
#include <cmath>
 
Implemente_instanciable(Op_Conv_EF_Stab_PolyMAC_HFV_Elem, "Op_Conv_EF_Stab_PolyMAC_HFV_Elem|Op_Conv_EF_Stab_PolyMAC_MPFA_Elem", Op_Conv_PolyMAC_CDO_base);
Implemente_instanciable_sans_constructeur(Op_Conv_Amont_PolyMAC_HFV_Elem, "Op_Conv_Amont_PolyMAC_HFV_Elem|Op_Conv_Amont_PolyMAC_MPFA_Elem", Op_Conv_EF_Stab_PolyMAC_HFV_Elem);
Implemente_instanciable_sans_constructeur(Op_Conv_Centre_PolyMAC_HFV_Elem, "Op_Conv_Centre_PolyMAC_HFV_Elem|Op_Conv_Centre_PolyMAC_MPFA_Elem", Op_Conv_EF_Stab_PolyMAC_HFV_Elem);
// XD Op_Conv_EF_Stab_PolyMAC_HFV_Elem interprete Op_Conv_EF_Stab_PolyMAC_HFV_Elem BRACE Class
// XD_CONT Op_Conv_EF_Stab_PolyMAC_HFV_Elem
 
Op_Conv_Amont_PolyMAC_HFV_Elem::Op_Conv_Amont_PolyMAC_HFV_Elem() { alpha_ = 1.0; }
Op_Conv_Centre_PolyMAC_HFV_Elem::Op_Conv_Centre_PolyMAC_HFV_Elem() { alpha_ = 0.0; }
 
Sortie& Op_Conv_EF_Stab_PolyMAC_HFV_Elem::printOn(Sortie& os) const { return Op_Conv_PolyMAC_CDO_base::printOn(os); }
Sortie& Op_Conv_Amont_PolyMAC_HFV_Elem::printOn(Sortie& os) const { return Op_Conv_EF_Stab_PolyMAC_HFV_Elem::printOn(os); }
Sortie& Op_Conv_Centre_PolyMAC_HFV_Elem::printOn(Sortie& os) const { return Op_Conv_EF_Stab_PolyMAC_HFV_Elem::printOn(os); }
 
Entree& Op_Conv_Amont_PolyMAC_HFV_Elem::readOn(Entree& is) { return Op_Conv_EF_Stab_PolyMAC_HFV_Elem::readOn(is); }
Entree& Op_Conv_Centre_PolyMAC_HFV_Elem::readOn(Entree& is) { return Op_Conv_EF_Stab_PolyMAC_HFV_Elem::readOn(is); }
 
Entree& Op_Conv_EF_Stab_PolyMAC_HFV_Elem::readOn(Entree& is)
{
  Op_Conv_PolyMAC_CDO_base::readOn(is);
  if (que_suis_je().debute_par("Op_Conv_EF_Stab") or que_suis_je().debute_par("Op_Conv_ALE")) //not in Op_Conv_Amont/Centre
    {
      Param param(que_suis_je());
      param.ajouter("alpha", &alpha_);            // XD_ADD_P double
      // XD_CONT parameter adjusting stabilization from 0 (centered scheme) to 1 (upwind scheme)
      param.lire_avec_accolades_depuis(is);
    }
 
  if (sub_type(Masse_Multiphase, equation())) //convection in Masse_Multiphase -> mass flow / title fields
    {
      const Pb_Multiphase& pb = ref_cast(Pb_Multiphase, equation().probleme());
      noms_cc_phases_.dimensionner(pb.nb_phases());
      cc_phases_.resize(pb.nb_phases());
 
      noms_vd_phases_.dimensionner(pb.nb_phases());
      vd_phases_.resize(pb.nb_phases());
 
      noms_x_phases_.dimensionner(pb.nb_phases());
      x_phases_.resize(pb.nb_phases());
 
      for (int i = 0; i < pb.nb_phases(); i++)
        {
          noms_cc_phases_[i] = Nom("debit_") + pb.nom_phase(i);
          noms_vd_phases_[i] = Nom("vitesse_debitante_") + pb.nom_phase(i);
          noms_x_phases_[i] = Nom("titre_") + pb.nom_phase(i);
        }
    }
  return is;
}
 
void Op_Conv_EF_Stab_PolyMAC_HFV_Elem::preparer_calcul()
{
  Op_Conv_PolyMAC_CDO_base::preparer_calcul();
 
  /* just in case... */
  const Domaine_Poly_base& domaine = le_dom_poly_.valeur();
  equation().init_champ_convecte();
  flux_bords_.resize(domaine.premiere_face_int(), (le_champ_inco ? le_champ_inco->valeurs() : equation().inconnue().valeurs()).line_size());
 
  if (domaine.domaine().nb_joints() && domaine.domaine().joint(0).epaisseur() < 2)
    {
      Cerr << "Op_Conv_EF_Stab_PolyMAC_HFV_Elem : largeur de joint insuffisante (minimum 2)!" << finl;
      Process::exit();
    }
}
 
double Op_Conv_EF_Stab_PolyMAC_HFV_Elem::calculer_dt_stab() const
{
  const DoubleTab& vit = vitesse_->valeurs();
  return calculer_dt_stab_gen(vit);
}
 
double Op_Conv_EF_Stab_PolyMAC_HFV_Elem::calculer_dt_stab_gen(const DoubleTab& vit) const
{
  const Domaine_Poly_base& domaine = le_dom_poly_.valeur();
  const DoubleVect& fs = domaine.face_surfaces(), &pf = equation().milieu().porosite_face(),
                    &ve = domaine.volumes(), &pe = equation().milieu().porosite_elem();
  const DoubleTab* alp = sub_type(Pb_Multiphase, equation().probleme()) ?
                         &ref_cast(Pb_Multiphase, equation().probleme()).equation_masse().inconnue().passe() : nullptr;
 
  const IntTab& e_f = domaine.elem_faces(), &f_e = domaine.face_voisins(), &fcl = ref_cast(Champ_Elem_PolyMAC_CDO, equation().inconnue()).fcl();
  const int N = std::min(vit.line_size(), equation().inconnue().valeurs().line_size());
 
  DoubleTrav flux(N); //sum of fluxes pf * |f| * vf
  double dt = 1.e10;
 
  for (int e = 0; e < domaine.nb_elem(); e++)
    {
      flux = 0.;
      const double poro_vol = pe(e) * ve(e);
      for (int i = 0; i < e_f.dimension(1); i++)
        {
          const int f = e_f(e, i);
          if (f < 0) continue;
 
          if (!Option_PolyMAC_family::TRAITEMENT_AXI || (Option_PolyMAC_family::TRAITEMENT_AXI && !(fcl(f, 0) == 4 || fcl(f, 0) == 5)))
            for (int n = 0; n < N; n++)
              {
                const int idx_phase = idx_phase_transportante_ > -1 ? idx_phase_transportante_ : n;
                flux(n) += pf(f) * fs(f) * std::max((e == f_e(f, 1) ? 1 : -1) * vit(f, idx_phase), 0.); //only the flux entering e counts
              }
        }
 
      for (int n = 0; n < N; n++)
        if ((!alp || (*alp)(e, n) > 1e-3) && std::abs(flux(n)) > 1e-12 /* avoid values 'very close to 0 but not completely zero' */)
          dt = std::min(dt, poro_vol / flux(n));
    }
 
  return Process::mp_min(dt);
}
 
void Op_Conv_EF_Stab_PolyMAC_HFV_Elem::dimensionner_blocs(matrices_t mats, const tabs_t& semi_impl) const
{
  const Domaine_Poly_base& domaine = le_dom_poly_.valeur();
  const Champ_Inc_base& cc = equation().champ_convecte();
 
  const IntTab& f_e = domaine.face_voisins(), &fcl_v = ref_cast(Champ_Face_base, vitesse_.valeur()).fcl();
  const int N = equation().inconnue().valeurs().line_size();
 
  for (const auto &i_m : mats)
    {
      // Check conditions on "vitesse" or semi-implicit derivatives
      if (i_m.first == "vitesse" || (cc.derivees().count(i_m.first) && !semi_impl.count(cc.le_nom().getString())))
        {
          Matrice_Morse mat;
          Stencil stencil(0, 2);
 
          int M = equation().probleme().get_champ(i_m.first.c_str()).valeurs().line_size();
 
          if (i_m.first == "vitesse")
            {
              // Specific treatment for "vitesse"
              for (int f = 0; f < domaine.nb_faces_tot(); f++)
                if (fcl_v(f, 0) < 2) // boundary condition but not Dirichlet, or internal face
                  for (int i = 0; i < 2; i++)
                    {
                      int e = f_e(f, i);
                      if (e < 0) continue;
 
                      if (e < domaine.nb_elem())
                        for (int n = 0; n < N; n++)
                          {
                            const int idx_phase = idx_phase_transportante_ > -1 ? idx_phase_transportante_ : n * (M > 1);
                            stencil.append_line(N * e + n, M * f + idx_phase);
                          }
                    }
            }
          else
            {
              // Treatment for scalar unknowns
              for (int f = 0; f < domaine.nb_faces_tot(); f++)
                for (int i = 0; i < 2; i++)
                  {
                    int e = f_e(f, i);
                    if (e < 0) continue;
 
                    if (e < domaine.nb_elem())
                      for (int j = 0; j < 2; j++)
                        {
                          int eb = f_e(f, j);
                          if (eb < 0) continue;
 
                          int m = 0;
                          for (int n = 0; n < N; n++, m += (M > 1))
                            stencil.append_line(N * e + n, M * eb + m);
                        }
                  }
            }
 
          // Remove duplicates and allocate the matrix
          tableau_trier_retirer_doublons(stencil);
          Matrix_tools::allocate_morse_matrix(equation().inconnue().valeurs().size_totale(),
                                              (i_m.first == "vitesse" ? vitesse_->valeurs().size_totale() : cc.derivees().at(i_m.first).size_totale()),
                                              stencil, mat);
 
          // Update the matrix in the collection
          if (i_m.second->nb_colonnes())
            *i_m.second += mat;
          else
            *i_m.second = mat;
        }
    }
}
 
void Op_Conv_EF_Stab_PolyMAC_HFV_Elem::ajouter_blocs(matrices_t mats, DoubleTab& secmem, const tabs_t& semi_impl) const
{
  const DoubleTab& vit = vitesse_->valeurs();
  ajouter_blocs_gen(mats, secmem, vit, semi_impl);
}
 
// adds the convection contribution to the right-hand side resu
// returns resu
void Op_Conv_EF_Stab_PolyMAC_HFV_Elem::ajouter_blocs_gen(matrices_t mats, DoubleTab& secmem, const DoubleTab& vit, const tabs_t& semi_impl) const
{
  const Domaine_Poly_base& domaine = le_dom_poly_.valeur();
  const Champ_Inc_base& cc = le_champ_inco ? le_champ_inco.valeur() : equation().champ_convecte();
  const std::string& nom_cc = cc.le_nom().getString();
  Matrice_Morse *m_vit = mats.count("vitesse") ? mats.at("vitesse") : nullptr;
 
  const IntTab& f_e = domaine.face_voisins(),
                &fcl = ref_cast(Champ_Inc_P0_base, equation().inconnue()).fcl(),
                 &fcl_v = ref_cast(Champ_Face_base, vitesse_.valeur()).fcl();
 
  const DoubleVect& fs = domaine.face_surfaces(), &pf = equation().milieu().porosite_face();
  const DoubleTab& vcc = semi_impl.count(nom_cc) ? semi_impl.at(nom_cc) : cc.valeurs(), bcc = cc.valeur_aux_bords();
  const int N = vcc.line_size(), Mv = vit.line_size();
 
  std::vector<std::tuple<const DoubleTab*, Matrice_Morse*, int>> d_cc; //list of cc derivatives to fill: pairs (cc derivative, matrix, number of variable components)
 
  if (!semi_impl.count(nom_cc))
    for (auto &i_m : mats)
      if (cc.derivees().count(i_m.first))
        d_cc.push_back(std::make_tuple(&cc.derivees().at(i_m.first), i_m.second, equation().probleme().get_champ(i_m.first.c_str()).valeurs().line_size()));
 
  DoubleTrav dv_flux(N), dc_flux(2, N); //derivatives of the convective flux at the face with respect to velocity / upwind / downwind convected field
 
  // Convection at internal faces, Neumann_val_ext or Dirichlet
  for (int f = 0; f < domaine.nb_faces(); f++)
    {
      const bool traiter_face = (fcl(f, 0) == 0 || (fcl(f, 0) > 4 && fcl(f, 0) < 7));
      if (traiter_face)
        {
          dv_flux = 0.;
          dc_flux = 0.;
 
          // Compute dv_flux and dc_flux
          for (int i = 0; i < 2; i++)
            {
              int e = f_e(f, i);
              int m = 0;
 
              for (int n = 0; n < N; n++, m += (Mv > 1))
                {
                  const int idx_phase = idx_phase_transportante_ > -1 ? idx_phase_transportante_ : m;
                  const double vit_f = vit(f, idx_phase);
                  const double v = vit_f ? vit_f : DBL_MIN;
                  const double fac = pf(f) * fs(f) * (1. + (v * (i ? -1 : 1) > 0 ? 1. : -1) * alpha_) / 2;
 
                  dv_flux(n) += fac * (e >= 0 ? vcc(e, n) : bcc(f, n));  // f is real -> trivial index in bcc
                  dc_flux(i, n) = e >= 0 ? fac * vit_f : 0;
                }
            }
 
          // Update the right-hand side
          for (int i = 0; i < 2; i++)
            {
              int e = f_e(f, i);
              if (e < 0) continue;
 
              if (e < domaine.nb_elem())
                {
                  int m = 0;
                  for (int n = 0; n < N; n++, m += (Mv > 1))
                    {
                      const int idx_phase = idx_phase_transportante_ > -1 ? idx_phase_transportante_ : m;
                      secmem(e, n) -= (i ? -1 : 1) * dv_flux(n) * vit(f, idx_phase);
                    }
                }
            }
 
          // Update derivatives: velocity
          if (m_vit && fcl_v(f, 0) < 2)
            for (int i = 0; i < 2; i++)
              {
                int e = f_e(f, i);
                if (e < 0) continue;
 
                if (e < domaine.nb_elem())
                  {
                    int m = 0;
                    for (int n = 0; n < N; n++, m += (Mv > 1))
                      {
                        const int idx_phase = idx_phase_transportante_ > -1 ? idx_phase_transportante_ : m;
                        (*m_vit)(N * e + n, Mv * f + idx_phase) += (i ? -1 : 1) * dv_flux(n);
                      }
                  }
              }
 
          // Update derivatives: convected field
          for (auto &&d_m_i : d_cc)
            {
              int M = std::get<2>(d_m_i);
              for (int i = 0; i < 2; i++)
                {
                  int e = f_e(f, i);
                  if (e < 0) continue;
 
                  if (e < domaine.nb_elem())
                    for (int j = 0; j < 2; j++)
                      {
                        int eb = f_e(f, j);
                        if (eb < 0) continue;
 
                        int m = 0;
                        for (int n = 0; n < N; n++, m += (M > 1))
                          (*std::get<1>(d_m_i))(N * e + n, M * eb + m) += (i ? -1 : 1) * dc_flux(j, n) * (*std::get<0>(d_m_i))(eb, m);
                      }
                }
            }
        }
    }
}
 
void Op_Conv_EF_Stab_PolyMAC_HFV_Elem::get_noms_champs_postraitables(Noms& nom, Option opt) const
{
  Op_Conv_PolyMAC_CDO_base::get_noms_champs_postraitables(nom, opt);
  Noms noms_compris;
 
  if (sub_type(Masse_Multiphase, equation()))
    {
      const Pb_Multiphase& pb = ref_cast(Pb_Multiphase, equation().probleme());
 
      for (int i = 0; i < pb.nb_phases(); i++)
        {
          noms_compris.add(noms_cc_phases_[i]);
          noms_compris.add(noms_vd_phases_[i]);
          noms_compris.add(noms_x_phases_[i]);
        }
    }
  if (opt == DESCRIPTION)
    Cerr << " Op_Conv_EF_Stab_PolyMAC_HFV_Elem : " << noms_compris << finl;
  else
    nom.add(noms_compris);
}
 
void Op_Conv_EF_Stab_PolyMAC_HFV_Elem::creer_champ(const Motcle& motlu)
{
  Op_Conv_PolyMAC_CDO_base::creer_champ(motlu);
 
  const int i = noms_cc_phases_.rang(motlu), j = noms_vd_phases_.rang(motlu), k = noms_x_phases_.rang(motlu);
 
  if (i >= 0 && !cc_phases_[i])
    {
      equation().discretisation().discretiser_champ("vitesse", equation().domaine_dis(), noms_cc_phases_[i], "kg/m2/s", dimension, 1, 0, cc_phases_[i]);
      champs_compris_.ajoute_champ(cc_phases_[i]);
    }
 
  if (j >= 0 && !vd_phases_[j])
    {
      equation().discretisation().discretiser_champ("vitesse", equation().domaine_dis(), noms_vd_phases_[j], "m/s", dimension, 1, 0, vd_phases_[j]);
      champs_compris_.ajoute_champ(vd_phases_[j]);
    }
 
  if (k >= 0 && !x_phases_[k])
    {
      equation().discretisation().discretiser_champ("temperature", equation().domaine_dis(), noms_x_phases_[k], "m/s", 1, 1, 0, x_phases_[k]);
      champs_compris_.ajoute_champ(x_phases_[k]);
    }
}
 
void Op_Conv_EF_Stab_PolyMAC_HFV_Elem::mettre_a_jour(double temps)
{
  const DoubleTab& vit = vitesse_->valeurs();
  mettre_a_jour_gen(temps, vit);
}
 
void Op_Conv_EF_Stab_PolyMAC_HFV_Elem::mettre_a_jour_gen(double temps, const DoubleTab& vit)
{
  Op_Conv_PolyMAC_CDO_base::mettre_a_jour(temps);
 
  const Domaine_Poly_base& domaine = le_dom_poly_.valeur();
  const Champ_Inc_base& cc = le_champ_inco ? le_champ_inco.valeur() : equation().champ_convecte();
 
  const IntTab& f_e = domaine.face_voisins(), &e_f = domaine.elem_faces();
 
  const DoubleVect& pf = equation().milieu().porosite_face(), &pe = equation().milieu().porosite_elem(),
                    &fs = domaine.face_surfaces(), &ve = domaine.volumes();
 
  const DoubleTab& vcc = cc.valeurs(),
                   &bcc = cc.valeur_aux_bords(), &xv = domaine.xv(), &xp = domaine.xp();
 
  DoubleTrav balp;
 
  if (vd_phases_.size())
    balp = equation().inconnue().valeur_aux_bords();
 
  const int  D = dimension,  N = vcc.line_size(), M = vit.line_size();
  DoubleTrav cc_f(N); //value of the convected field at faces
 
  /* boundary fluxes */
  for (int f = 0; f < domaine.premiere_face_int(); f++)
    {
      cc_f = 0.;
      for (int i = 0; i < 2; i++)
        {
          int e = f_e(f, i);
          int m = 0;
          for (int n = 0; n < N; n++, m += (M > 1))
            {
              const int idx_phase = idx_phase_transportante_ > -1 ? idx_phase_transportante_ : m;
              cc_f(n) += (1. + (vit(f, idx_phase) * (i ? -1 : 1) >= 0 ? 1. : -1.) * alpha_) / 2 * (e >= 0 ? vcc(e, n) : bcc(f, n));
            }
        }
 
      int m = 0;
      for (int n = 0; n < N; n++, m += (M > 1))
        {
          const int idx_phase = idx_phase_transportante_ > -1 ? idx_phase_transportante_ : m;
          flux_bords_(f, n) = pf(f) * fs(f) * vit(f, idx_phase) * cc_f(n);
        }
    }
 
  if (cc_phases_.size())
    {
      int m = 0;
      for (int n = 0; n < N; n++, m += (M > 1))
        if (cc_phases_[n]) /* update flow rate fields */
          {
            Champ_Face_PolyMAC_HFV& c_ph = ref_cast(Champ_Face_PolyMAC_HFV, cc_phases_[n].valeur());
            DoubleTab& v_ph = c_ph.valeurs();
 
            for (int f = 0; f < domaine.nb_faces(); f++)
              {
                v_ph(f) = 0.;
                const int idx_phase = idx_phase_transportante_ > -1 ? idx_phase_transportante_ : m;
                if (m != idx_phase) continue;
                for (int i = 0; i < 2; i++)
                  {
                    int e = f_e(f, i);
                    double signe = (vit(f, m) * (i ? -1 : 1) >= 0) ? 1. : -1.;
                    double facteur = (1. + signe * alpha_) / 2.;
 
                    v_ph(f) += facteur * (e >= 0 ? vcc(e, n) : bcc(f, n));
                  }
 
                v_ph(f) *= vit(f, m) * pf(f);
              }
 
            c_ph.changer_temps(temps);
          }
    }
 
  if (vd_phases_.size())
    {
      int m = 0;
      for (int n = 0; n < N; n++, m += (M > 1))
        if (vd_phases_[n]) /* update superficial velocity fields */
          {
            const DoubleTab& alp = equation().inconnue().valeurs();
            Champ_Face_PolyMAC_HFV& c_ph = ref_cast(Champ_Face_PolyMAC_HFV, vd_phases_[n].valeur());
            DoubleTab& v_ph = c_ph.valeurs();
 
            /* fill the face part, then ask the field to interpolate to elements */
            for (int f = 0; f < domaine.nb_faces(); f++)
              {
                v_ph(f) = 0.;
                const int idx_phase = idx_phase_transportante_ > -1 ? idx_phase_transportante_ : m;
                if (m != idx_phase) continue;
                for (int i = 0; i < 2; i++)
                  {
                    int e = f_e(f, i);
                    double signe = (vit(f, m) * (i ? -1 : 1) >= 0) ? 1. : -1.;
                    double facteur = (1. + signe * alpha_) / 2.;
 
                    v_ph(f) += facteur * (e >= 0 ? alp(e, n) : balp(f, n));
 
                  }
                v_ph(f) *= vit(f, m) * pf(f);
              }
 
            c_ph.changer_temps(temps);
          }
    }
 
  DoubleTrav G(N), v(N, D);
  double Gt;
 
  if (x_phases_.size())
    for (int e = 0; e < domaine.nb_elem(); e++) // void fraction: at elements
      {
        v = 0.;
        for (int i = 0; i < e_f.dimension(1); i++)
          {
            int f = e_f(e, i);
            if (f < 0) continue;
 
            for (int n = 0; n < N; n++)
              for (int d = 0; d < D; d++)
                v(n, d) += fs(f) * pf(f) * (xv(f, d) - xp(e, d)) * (e == f_e(f, 0) ? 1 : -1) * vit(f, n) / (pe(e) * ve(e));
          }
 
        Gt = 0.;
        for (int n = 0; n < N; n++)
          {
            G(n) = vcc(e, n) * sqrt(domaine.dot(&v(n, 0), &v(n, 0)));
            Gt += G(n);
          }
 
        for (int n = 0; n < N; n++)
          if (x_phases_[n])
            x_phases_[n]->valeurs()(e) = Gt ? G(n) / Gt : 0.;
      }
 
  if (x_phases_.size())
    for (int n = 0; n < N; n++)
      if (x_phases_[n])
        x_phases_[n]->changer_temps(temps);
}