Op_Conv_VEF_Face.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_VEF_Face.h>
#include <Champ_P1NC.h>
#include <Porosites_champ.h>
#include <Convection_tools.h>
#include <Dirichlet_homogene.h>
#include <Periodique.h>
#include <Symetrie.h>
#include <Neumann_sortie_libre.h>
#include <TRUSTVect.h>
#include <Device.h>
#include <Tetra_VEF.h>
#include <Tri_VEF.h>
#include <Comm_Group_MPI.h>
 
 
 
//If constexpr can be very useful for performances (see https://rbourgeois33.github.io./posts/post1/#template-away-heavy-branches)
//But it is not compatible with old compilers. We fallback to plain if when necessary
#if defined(__cpp_if_constexpr)
#define TRUST_IFCONSTEXPR if constexpr
#else
#define TRUST_IFCONSTEXPR if
#pragma message("Warning: your C++ std is old and does not enables 'if constexpr', switching to plain 'if'. Think about upgrading to C++ 17.")
#endif
 
Implemente_instanciable_sans_constructeur(Op_Conv_VEF_Face,"Op_Conv_Generic_VEF_P1NC",Op_Conv_VEF_base);
// XD convection_generic convection_deriv generic NO_BRACE Keyword for generic calling of upwind and muscl convective
// XD_CONT scheme in VEF discretization. For muscl scheme, limiters and order for fluxes calculations have to be
// XD_CONT specified. The available limiters are : minmod - vanleer -vanalbada - chakravarthy - superbee, and the order
// XD_CONT of accuracy is 1 or 2. Note that chakravarthy is a non-symmetric limiter and superbee may engender results
// XD_CONT out of physical limits. By consequence, these two limiters are not recommended. NL2 Examples: NL2 convection
// XD_CONT { generic amont }NL2 convection { generic muscl minmod 1 }NL2 convection { generic muscl vanleer 2 }NL2 NL2
// XD_CONT In case of results out of physical limits with muscl scheme (due for instance to strong non-conformal
// XD_CONT velocity flow field), user can redefine in data file a lower order and a smoother limiter, as : convection {
// XD_CONT generic muscl minmod 1 }
// XD attr type chaine(into=["amont","muscl","centre"]) type REQ type of scheme
// XD attr limiteur chaine(into=["minmod","vanleer","vanalbada","chakravarthy","superbee"]) limiteur OPT type of limiter
// XD attr ordre entier(into=[1,2,3]) ordre OPT order of accuracy
// XD attr alpha floattant alpha OPT alpha
 
Sortie& Op_Conv_VEF_Face::printOn(Sortie& s ) const
{
  return s << que_suis_je() ;
}
 
Entree& Op_Conv_VEF_Face::readOn(Entree& s )
{
  Motcle type_op_lu;
  s >> type_op_lu;
 
  if (!(type_op_lu=="amont") && !(type_op_lu=="muscl") && !(type_op_lu=="centre"))
    {
      Cerr << type_op_lu << " is not understood by " << que_suis_je() << finl;
      Cerr << " choose among: amont - muscl - centre " << finl;
      exit();
    }
 
  if (type_op_lu=="muscl")
    {
      type_op = muscl;
      s >> type_lim;
 
      if ( !(type_lim=="minmod") && !(type_lim=="vanleer") && !(type_lim=="vanalbada")
           &&  !(type_lim=="chakravarthy") && !(type_lim=="superbee") )
        {
          Cerr << type_lim << " is not understood by " << que_suis_je() << finl;
          Cerr << " choose among: minmod - vanleer - vanalbada - chakravarthy - superbee " << finl;
          exit();
        }
 
      if (type_lim=="minmod")
        {
          LIMITEUR=&minmod;
          type_lim_int = type_lim_minmod;
        }
      if (type_lim=="vanleer")
        {
          LIMITEUR=&vanleer;
          type_lim_int = type_lim_vanleer;
        }
      if (type_lim=="vanalbada")
        {
          LIMITEUR=&vanalbada;
          type_lim_int = type_lim_vanalbada;
        }
      if (type_lim=="chakravarthy")
        {
          LIMITEUR=&chakravarthy;
          type_lim_int = type_lim_chakravarthy;
        }
      if (type_lim=="superbee")
        {
          LIMITEUR=&superbee;
          type_lim_int = type_lim_superbee;
        }
      s >> ordre_;
 
      if (ordre_!=1 && ordre_!=2 && ordre_!=3)
        {
          Cerr << "the order after " << type_lim << " in " << que_suis_je() << " must be 1, 2, or 3" <<  finl;
          exit();
        }
      if (ordre_==3)
        {
          // Read alpha_
          s >> alpha_;
        }
    }
 
  if (type_op_lu=="centre")
    {
      type_op = centre;
      s >> ordre_;
 
      if (ordre_!=1 && ordre_!=2)
        {
          Cerr << "the order after " << type_op_lu << " in " << que_suis_je() << " must be 1 or 2 " <<  finl;
          exit();
        }
    }
 
  if (type_op_lu=="amont")
    {
      type_op = amont;
      ordre_ = 1;
    }
 
  return s ;
}
 
void Op_Conv_VEF_Face::completer()
{
  Op_Conv_VEF_base::completer();
  // For an element that is not on the boundary (rang==-1),
  // check if one of its vertices is on the boundary (vertex array) (C MALOD 17/07/2007)
  if (type_elem_Cl_.size_array() == 0)
    {
      const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, le_dom_vef.valeur());
      const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
      const int nb_elem_tot = domaine_VEF.nb_elem_tot();
      const IntVect& rang_elem_non_std = domaine_VEF.rang_elem_non_std();
      type_elem_Cl_.resize(nb_elem_tot);
      for (int poly = 0; poly < nb_elem_tot; poly++)
        {
          int rang = rang_elem_non_std(poly);
          type_elem_Cl_[poly] = rang == -1 ? 0 : domaine_Cl_VEF.type_elem_Cl(rang);
        }
    }
}
 
struct FluxTetraKernelData
{
  // Scalars
  int nb_elem;
  int nb_elem_tot;
  int nb_faces;
  int nb_faces_bord;
  int nfa7;
  int ncomp_ch_transporte;
  double alpha;
  int marq;
  bool option_calcul_flux_en_un_point;
  bool option_appliquer_cl_dirichlet;
  bool isMuscl;
  bool isAmont;
 
  // Views
  CIntArrView rang_elem_non_std_v;
  CIntTabView elem_faces_v;
  CDoubleArrView porosite_face_v;
  CDoubleArrView porosite_elem_v;
  CDoubleTabView coord_sommets_v;
  CIntTabView les_elems_v;
  CDoubleTabView3 facette_normales_v;
  CIntArrView est_une_face_de_dirichlet_v;
  CDoubleTabView xp_v;
  CDoubleTabView xv_v;
  CIntArrView type_elem_Cl_v;
  CIntArrView traitement_pres_bord_v;
  CIntTabView KEL_v;
  CDoubleTabView3 normales_facettes_Cl_v;
  CDoubleTabView4 vecteur_face_facette_Cl_v;
  CDoubleTabView vitesse_v;
  CDoubleTabView transporte_face_v;
  CDoubleTabView3 gradient_v;
 
  // Outputs
  DoubleTabView resu_v;
  DoubleTabView flux_b_v;
};
 
template<int ordre, bool isMuscl, bool virtual_only>
void compute_flux_tetra_kernel(const FluxTetraKernelData& kernel_data)
{
 
  //const int dim = Objet_U::dimension;
  const int dim = 3; // Help compiler: +16% speed-up ! Not the same effect if setting nsom=4 and nfac=4 for a Tetra
  const int nfac_ = 4;
  const int nsom_ = 4;
  const double third = 1.0/3;
  const double twelvth = 1.0/(3+3*3);
 
  // Unpack scalars
  int nb_faces                     = kernel_data.nb_faces;
  int nb_faces_bord                = kernel_data.nb_faces_bord;
  //int nfa7                         = kernel_data.nfa7;
  int ncomp_ch_transporte          = kernel_data.ncomp_ch_transporte;
  double alpha                     = kernel_data.alpha;
  int marq                          = kernel_data.marq;
  bool option_calcul_flux_en_un_point     = kernel_data.option_calcul_flux_en_un_point;
  bool option_appliquer_cl_dirichlet      = kernel_data.option_appliquer_cl_dirichlet;
  bool isAmont                     = kernel_data.isAmont;
 
  // Unpack views
  CIntArrView rang_elem_non_std_v           = kernel_data.rang_elem_non_std_v;
  CIntTabView elem_faces_v                  = kernel_data.elem_faces_v;
  CDoubleArrView porosite_face_v            = kernel_data.porosite_face_v;
  CDoubleArrView porosite_elem_v            = kernel_data.porosite_elem_v;
  CDoubleTabView coord_sommets_v           = kernel_data.coord_sommets_v;
  CIntTabView les_elems_v                   = kernel_data.les_elems_v;
  CDoubleTabView3 facette_normales_v       = kernel_data.facette_normales_v;
  CIntArrView est_une_face_de_dirichlet_v  = kernel_data.est_une_face_de_dirichlet_v;
  CDoubleTabView xp_v                       = kernel_data.xp_v;
  CDoubleTabView xv_v                       = kernel_data.xv_v;
  CIntArrView type_elem_Cl_v                = kernel_data.type_elem_Cl_v;
  CIntArrView traitement_pres_bord_v       = kernel_data.traitement_pres_bord_v;
  CIntTabView KEL_v                         = kernel_data.KEL_v;
  CDoubleTabView3 normales_facettes_Cl_v   = kernel_data.normales_facettes_Cl_v;
  CDoubleTabView4 vecteur_face_facette_Cl_v= kernel_data.vecteur_face_facette_Cl_v;
  CDoubleTabView vitesse_v                  = kernel_data.vitesse_v;
  CDoubleTabView transporte_face_v         = kernel_data.transporte_face_v;
  CDoubleTabView3 gradient_v               = kernel_data.gradient_v;
 
  // Unpack outputs
  DoubleTabView resu_v                     = kernel_data.resu_v;
  DoubleTabView flux_b_v                   = kernel_data.flux_b_v;
 
  //Only virtual elements
  const int nbeg = virtual_only ? kernel_data.nb_elem : 0;
  const int nend = virtual_only ? kernel_data.nb_elem_tot : kernel_data.nb_elem;
 
  //The reel kernel does not start the gpu counter because it is manager by the echange_espace_virtuel_async
  //Only the virtual one does because it is synchronous, after the comms
  if (virtual_only)
    {
      start_gpu_timer(__KERNEL_NAME__);
    }
 
#ifdef TRUST_USE_GPU
  Kokkos::parallel_for(Kokkos::MDRangePolicy<Kokkos::Rank<2>>({nbeg,0}, {nend, kernel_data.nfa7}), KOKKOS_LAMBDA(const int poly, const int fa7)
#else
  Kokkos::parallel_for(Kokkos::RangePolicy<>(nbeg, nend), KOKKOS_LAMBDA(const int poly)
#endif
  {
    int contrib = 0;
    // compute the face indices of the polyhedron
    int face[4];
    for (int face_adj = 0; face_adj < nfac_; face_adj++)
      {
        int fac = elem_faces_v(poly, face_adj);
        face[face_adj] = fac;
        if (fac < nb_faces) contrib = 1; // A real face on the virtual element
      }
    //
    if (contrib)
      {
        //prefetch as much as possible
        bool calcul_flux_en_un_point = (ordre != 3) && (ordre == 1 || traitement_pres_bord_v(poly));
        bool flux_en_un_point = calcul_flux_en_un_point || option_calcul_flux_en_un_point;
        double vs[3] = {0.0, 0.0, 0.0};
        double vsom[12];  // 4 nodes × 3 components
        double poro[4];
        double vfa[12];
        int les_elems_[4] = { les_elems_v(poly,0), les_elems_v(poly,1), les_elems_v(poly,2), les_elems_v(poly,3) };
        // Porosity handling
        double coeff = (marq == 0) ? 1. / porosite_elem_v(poly) : 1.0;
        int itypcl = type_elem_Cl_v(poly);
        // Determine the BC type based on the rank
        int rang = rang_elem_non_std_v(poly);
        double xc[3];
        TRUST_IFCONSTEXPR (ordre == 3) // To be optimised! Risk of wrong results in parallel if ordre=3
        {
          double xsom[12];
          for (int i = 0; i < nsom_; i++)
            for (int j = 0; j < dim; j++)
              xsom[i * 3 + j] = coord_sommets_v(les_elems_[i], j);
          int idirichlet, n1, n2, n3;
          calcul_xg_tetra(xc, xsom, itypcl, idirichlet, n1, n2, n3);
        }
 
        double xp[3] = { xp_v(poly,0), xp_v(poly,1), xp_v(poly,2) };
 
        // Porosity + velocity contributions
        for (int i = 0; i < nfac_; i++)
          {
            poro[i] = porosite_face_v(face[i]);
 
            for (int j = 0; j < dim; j++)
              {
                double vf = vitesse_v(face[i], j);
                vfa[i*3 + j] = vf;
                vs[j] += vf * poro[i];
              }
          }
        for (int i = 0; i < nfac_; i++)
          {
            for (int j = 0; j < dim; j++)
              {
                vsom[i*3 + j] = vs[j] - 3.0 * vfa[i*3 + j] * poro[i];
              }
          }
 
        double vc[3];
        calcul_vc_tetra(face, vc, vs, vsom, vfa, itypcl, poro);
        // compute xc (at the intersection of the 3 facets), needed for muscl3
 
#ifndef TRUST_USE_GPU
        for (int fa7 = 0; fa7 < kernel_data.nfa7; fa7++)
#endif
          {
            int KEL_v_[4] = {KEL_v(0, fa7), KEL_v(1, fa7),KEL_v(2, fa7),KEL_v(3, fa7)};
            int les_kel_elems_[2]= {0,0};
 
            int num10 = -1;
            int num20 = -1;
 
            double psc_c = 0, psc_s = 0, psc_s2 = 0;
 
            // facet normals
            double cc[3];
            for (int j = 0; j < 3; j++)
              {
                cc[j] = (rang == -1)
                        ? facette_normales_v(poly, fa7, j)
                        : normales_facettes_Cl_v(rang, fa7, j);
              }
 
            // Compute velocities at C,S,S2 (the 3 endpoints of fa7) and M (the center of fa7)
            for (int i = 0; i < 3; i++)
              {
                psc_c += vc[i] * cc[i];
              }
 
            switch (KEL_v_[0])
              {
              case 0:
                num10 = face[0];
                break;
              case 1:
                num10 = face[1];
                break;
              case 2:
                num10 = face[2];
                break;
              case 3:
                num10 = face[3];
 
                break;
              default: /* handle error if needed */
                break;
              }
            switch (KEL_v_[1])
              {
              case 0:
                num20 = face[0];
                break;
              case 1:
                num20 = face[1];
                break;
              case 2:
                num20 = face[2];
                break;
              case 3:
                num20 = face[3];
                break;
              default: /* handle error if needed */
                break;
              }
            switch (KEL_v_[2])
              {
              case 0:
                psc_s += vsom[0] * cc[0] + vsom[1] * cc[1] + vsom[2] * cc[2];
                les_kel_elems_[0]=les_elems_[0];
 
                break;
              case 1:
                psc_s += vsom[3] * cc[0] + vsom[4] * cc[1] + vsom[5] * cc[2];
                les_kel_elems_[0]=les_elems_[1];
                break;
              case 2:
                psc_s += vsom[6] * cc[0] + vsom[7] * cc[1] + vsom[8] * cc[2];
                les_kel_elems_[0]=les_elems_[2];
 
                break;
              case 3:
                psc_s += vsom[9] * cc[0] + vsom[10] * cc[1] + vsom[11] * cc[2];
                les_kel_elems_[0]=les_elems_[3];
 
                break;
              }
            switch (KEL_v_[3])
              {
              case 0:
                psc_s2 += vsom[0] * cc[0] + vsom[1] * cc[1] + vsom[2] * cc[2];
                les_kel_elems_[1]=les_elems_[0];
                break;
              case 1:
                psc_s2 += vsom[3] * cc[0] + vsom[4] * cc[1] + vsom[5] * cc[2];
                les_kel_elems_[1]=les_elems_[1];
                break;
              case 2:
                psc_s2 += vsom[6] * cc[0] + vsom[7] * cc[1] + vsom[8] * cc[2];
                les_kel_elems_[1]=les_elems_[2];
                break;
              case 3:
                psc_s2 += vsom[9] * cc[0] + vsom[10] * cc[1] + vsom[11] * cc[2];
                les_kel_elems_[1]=les_elems_[3];
                break;
              }
 
            int sommet_s = les_kel_elems_[0];
            int sommet_s2 =  les_kel_elems_[1];
            psc_c *=coeff;
            psc_s *=coeff;
            psc_s2 *=coeff;
            double psc_m = (psc_c + psc_s + psc_s2) *third;
            // Apply Dirichlet BCs if num1 or num2 is a face with a Dirichlet BC:
            // in that case fa7 coincides with face num1 or num2 -> C is at the center of the face
 
            int appliquer_cl_dirichlet = 0;
            if (option_appliquer_cl_dirichlet)
              if (est_une_face_de_dirichlet_v(num10) || est_une_face_de_dirichlet_v(num20))
                {
                  appliquer_cl_dirichlet = 1;
                  psc_m = psc_c;
                }
 
            // Determine the upwind faces for points M,C,S,S2
            int face_amont_m = (psc_m >= 0) ? num10 : num20;
 
            [[maybe_unused]] int face_amont_c;
            int face_amont_s, face_amont_s2;
            TRUST_IFCONSTEXPR (ordre == 3 && isMuscl)
            {
              face_amont_c = ((psc_c >= 0) ? num10 : num20);
              face_amont_s = ((psc_s >= 0) ? num10 : num20) ;
              face_amont_s2 =  ((psc_s2 >= 0) ? num10 : num20) ;
            }
            else
              {
                face_amont_c= face_amont_m;
                face_amont_s= face_amont_m;
                face_amont_s2= face_amont_m;
              }
 
            int item_m, item_s, item_s2;
            [[maybe_unused]] int item_c;
 
            TRUST_IFCONSTEXPR (isMuscl)
            {
              // Use the "face amont" scheme (Muscl)
              item_m  = face_amont_m;
              item_c  = face_amont_c;
              item_s  = face_amont_s;
              item_s2 = face_amont_s2;
            }
            else
              {
                // Use the center-of-element scheme (non-Muscl)
                item_m  = poly;
                item_c  = poly;
                item_s  = poly;
                item_s2 = poly;
              }
 
            int dir =  (psc_m >= 0) ? 0 : 1;
            int kel =  (psc_m >= 0) ? KEL_v_[0] : KEL_v_[1];
            int num_face=0;
            switch (kel)
              {
              case 0:
                num_face = face[0];
                break;
              case 1:
                num_face = face[1];
                break;
              case 2:
                num_face = face[2];
                break;
              case 3:
                num_face = face[3];
                break;
              }
 
            double centre_fa7[3];
            if (rang == -1)
              {
                int isom_glob_0 = les_kel_elems_[0];//2
                int isom_glob_1 = les_kel_elems_[1];//3
 
                for (int j = 0; j < 3; j++)
                  {
                    centre_fa7[j] = xp[j]
                                    + coord_sommets_v(isom_glob_0, j)
                                    + coord_sommets_v(isom_glob_1, j);
                    centre_fa7[j] *= third;
                  }
              }
            bool flux_avec_m = isAmont || appliquer_cl_dirichlet;
 
            double xv_m_arr[3];
            double xv_s_arr[3];
            double xv_s2_arr[3];
            double cs_s_arr[3];
            double cs_s2_arr[3];
 
            for (int j = 0; j < 3; j++)
              {
                xv_m_arr[j]  = xv_v(num_face, j);
                xv_s_arr[j]  = xv_v(face_amont_s, j);
                xv_s2_arr[j] = xv_v(face_amont_s2, j);
                cs_s_arr[j]  = coord_sommets_v(sommet_s, j);
                cs_s2_arr[j] = coord_sommets_v(sommet_s2, j);
              }
 
            for (int comp0 = 0; comp0 < ncomp_ch_transporte; comp0++)
              {
                double inco_m = transporte_face_v(face_amont_m, comp0);
                double flux;
                if (flux_avec_m) // amont
                  flux = inco_m * psc_m;
                else // muscl or centre
                  {
                    // PL: data locality matters ! Try to access arrays one by one and store into registers
                    double inco_s = transporte_face_v(face_amont_s, comp0);
                    double inco_s2 =  transporte_face_v(face_amont_s2, comp0) ;
                    for (int j = 0; j < dim; j++)
                      {
                        double gm = gradient_v(item_m, comp0, j);
                        double gs = gradient_v(item_s, comp0, j);
                        double gs2 =  gradient_v(item_s2, comp0, j);
 
                        // Compute the unknown at the center M of fa7
 
                        inco_m += gm  * (rang == -1 ? centre_fa7[j] - xv_m_arr[j]
                                         : vecteur_face_facette_Cl_v(rang, fa7, j, dir));
                        inco_s += gs  * (cs_s_arr[j]  - xv_s_arr[j]);
                        inco_s2+= gs2 * (cs_s2_arr[j] - xv_s2_arr[j]);
                      }
                    // Compute the unknown at C, the other extremity of fa7, intersection with the other fa7s
                    // of the polyhedron. C=G center of the polyhedron if extended volume is not used
                    // xc given by elemvef.calcul_xg()
                    double inco_c;
                    TRUST_IFCONSTEXPR (ordre == 3)
                    {
                      inco_c = transporte_face_v(face_amont_c, comp0);
                      for (int j = 0; j < dim; j++)
                        inco_c += gradient_v(item_c, comp0, j) * (-xv_v(face_amont_c, j) + xc[j]);
                    }
                    else
                      inco_c = dim * inco_m - inco_s - inco_s2;
 
                    if (flux_en_un_point)
                      // Compute flux at 1 point
                      flux = inco_m * psc_m;
                    else
                      // Compute flux at 3 points
                      flux = (inco_c * psc_c + inco_s * psc_s + inco_s2 * psc_s2 + 9 * inco_m * psc_m)*twelvth;
                  }
 
                // Weighting by coefficient alpha
                flux *= alpha;
 
                int compo = ncomp_ch_transporte == 1 ? 0 : comp0;
                Kokkos::atomic_sub(&resu_v(num10, compo), flux);
                Kokkos::atomic_add(&resu_v(num20, compo), flux);
                if (num10 < nb_faces_bord)
                  Kokkos::atomic_add(&flux_b_v(num10, compo), flux);
                if (num20 < nb_faces_bord)
                  Kokkos::atomic_sub(&flux_b_v(num20, compo), flux);
 
              }// loop over components
          } // end of loop over facets
      } // end of loop
  });
  if (virtual_only) end_gpu_timer(__KERNEL_NAME__);
}
 
//Wrapper to hide the async logic
//Same name as the actual kernel function to appear as the same on the profiler
template<int ordre, bool isMuscl>
void compute_flux_tetra_kernel(const FluxTetraKernelData& kernel_data, DoubleTab& tab_gradient)
{
  //Packing, isend, irecv but no wait and no unpacking
  //If GPU, the isend/irdc are cudaMemcpyAsync P2P, on a given stream, no blocking !
  tab_gradient.start_echange_espace_virtuel_async(__KERNEL_NAME__);
  //Only operate on reel elements
  compute_flux_tetra_kernel<ordre, /*muscl*/ isMuscl, /* virtual poly only */ false>(kernel_data);
  //Wait and unpack virtual elements
  tab_gradient.finish_echange_espace_virtuel_async(__KERNEL_NAME__);
  //Only operates on virtual elements
  compute_flux_tetra_kernel<ordre, /*muscl*/ isMuscl, /* virtual poly only */ true>(kernel_data);
}
//
//   Methods of the class Op_Conv_VEF_Face
//
////////////////////////////////////////////////////////////////////
//
//                      Implementation of methods
//
//                   of the class Op_Conv_VEF_Face
//
////////////////////////////////////////////////////////////////////
 
DoubleTab& Op_Conv_VEF_Face::ajouter(const DoubleTab& transporte, DoubleTab& resu) const
{
  const Champ_Inc_base& la_vitesse_aux_faces = vitesse();
  return ajouter_gen(transporte, la_vitesse_aux_faces, resu);
}
 
DoubleTab& Op_Conv_VEF_Face::ajouter_gen(const DoubleTab& transporte, const Champ_Inc_base& la_vitesse, DoubleTab& resu) const
{
 
  assert((type_op==amont) || (type_op==muscl) || (type_op==centre));
  const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
  const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, le_dom_vef.valeur());
  const DoubleTab& velocity_tab=la_vitesse.valeurs();
  const DoubleVect& porosite_face = equation().milieu().porosite_face();
 
  int marq=phi_u_transportant(equation());
  DoubleTab transporte_face_;
  DoubleTab vitesse_face_;
  // either transporte_face=phi*transporte and vitesse_face=vitesse
  // or transporte_face=transporte and vitesse_face=phi*vitesse
  // depending on whether transport is done with phi*u or with u.
  const DoubleTab& transporte_face = modif_par_porosite_si_flag(transporte,transporte_face_,!marq,porosite_face);
  const DoubleTab& vitesse_face    = modif_par_porosite_si_flag(la_vitesse.valeurs(),vitesse_face_,marq,porosite_face);
 
  const IntTab& elem_faces = domaine_VEF.elem_faces();
  const auto& facette_normales = domaine_VEF.facette_normales();
  const Domaine& domaine = domaine_VEF.domaine();
  const int nfa7 = domaine_VEF.type_elem().nb_facette();
  const int nb_elem_tot = domaine_VEF.nb_elem_tot();
  const int nb_elem = domaine_VEF.nb_elem();
  const IntVect& rang_elem_non_std = domaine_VEF.rang_elem_non_std();
  const DoubleTab& normales_facettes_Cl = domaine_Cl_VEF.normales_facettes_Cl();
  int premiere_face_int = domaine_VEF.premiere_face_int();
  int nb_faces = domaine_VEF.nb_faces();
  int nfac = domaine.nb_faces_elem();
  int nsom = domaine.nb_som_elem();
  const DoubleTab& vecteur_face_facette_Cl = domaine_Cl_VEF.vecteur_face_facette_Cl();
  int nb_bord = domaine_VEF.nb_front_Cl();
  const IntTab& les_elems=domaine.les_elems();
 
  // Ensures a flux_bord consistent with the BCs (but sometimes diverges?)
  // Active only for order 3
  int option_appliquer_cl_dirichlet = (ordre_==3 ? 1 : 0);
  int option_calcul_flux_en_un_point = 0;//(ordre==3 ? 1 : 0);
  // Define an array for special treatment of schemes near boundaries
  if (traitement_pres_bord_.size_array()!=nb_elem_tot)
    {
      traitement_pres_bord_.resize_array(nb_elem_tot);
      traitement_pres_bord_=0;
      // For muscl3, apply the minmod limiter on elements with a Dirichlet face
      if (ordre_==3)
        {
          for (int n_bord=0; n_bord<nb_bord; n_bord++)
            {
              const Cond_lim_base& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord).valeur();
              if ( sub_type(Dirichlet,la_cl) || sub_type(Dirichlet_homogene,la_cl) )
                {
                  const Front_VF& le_bord = ref_cast(Front_VF,la_cl.frontiere_dis());
                  int nb_faces_tot = le_bord.nb_faces_tot();
                  const IntTab& face_voisins = domaine_VEF.face_voisins();
                  for (int ind_face=0; ind_face<nb_faces_tot; ind_face++)
                    {
                      int num_face = le_bord.num_face(ind_face);
                      int elem = face_voisins(num_face,0);
                      traitement_pres_bord_[elem]=1;
                    }
                }
            }
        }
      else
        {
          // For the current muscl/centre schemes, a first-order flux computation is used
          // at boundary cells or at cells having a Dirichlet vertex
          ArrOfInt est_un_sommet_de_bord_(domaine_VEF.nb_som_tot());
          for (int n_bord=0; n_bord<nb_bord; n_bord++)
            {
              const Cond_lim_base& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord).valeur();
              if ( sub_type(Dirichlet,la_cl) || sub_type(Dirichlet_homogene,la_cl) )
                {
                  const Front_VF& le_bord = ref_cast(Front_VF,la_cl.frontiere_dis());
                  int nb_faces_tot = le_bord.nb_faces_tot();
                  int size = domaine_VEF.face_sommets().dimension(1);
                  for (int ind_face=0; ind_face<nb_faces_tot; ind_face++)
                    for (int som=0; som<size; som++)
                      {
                        int face = le_bord.num_face(ind_face);
                        est_un_sommet_de_bord_[domaine_VEF.face_sommets(face,som)]=1;
                      }
                }
            }
          for (int elem=0; elem<nb_elem_tot; elem++)
            {
              if (rang_elem_non_std(elem)!=-1)
                traitement_pres_bord_[elem]=1;
              else
                {
                  for (int n_som=0; n_som<nsom; n_som++)
                    if (est_un_sommet_de_bord_[les_elems(elem,n_som)])
                      traitement_pres_bord_[elem]=1;
                }
            }
        }
      // Build the array est_une_face_de_dirichlet_
      est_une_face_de_dirichlet_.resize_array(domaine_VEF.nb_faces_tot());
      est_une_face_de_dirichlet_=0;
      for (int n_bord=0; n_bord<nb_bord; n_bord++)
        {
          const Cond_lim_base& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord).valeur();
          if ( sub_type(Dirichlet,la_cl) || sub_type(Dirichlet_homogene,la_cl) )
            {
              const Front_VF& le_bord = ref_cast(Front_VF,la_cl.frontiere_dis());
              int nb_faces_tot = le_bord.nb_faces_tot();
              for (int ind_face=0; ind_face<nb_faces_tot; ind_face++)
                {
                  int num_face = le_bord.num_face(ind_face);
                  est_une_face_de_dirichlet_[num_face] = 1;
                }
            }
        }
    }
 
  // For convection treatment, we distinguish standard polyhedra
  // which do not "see" boundary conditions, from non-standard
  // polyhedra which have at least one boundary face.
  // A standard polyhedron has n facets on which the convection
  // scheme is applied.
  // For a non-standard polyhedron with Dirichlet boundary conditions,
  // some facets are carried by boundary faces.
  // In short, for a polyhedron, convection treatment depends
  // on the type (triangle, tetrahedron ...) and the number of Dirichlet faces.
 
  const Elem_VEF_base& type_elemvef= domaine_VEF.type_elem();
  int istetra=0;
  Nom nom_elem=type_elemvef.que_suis_je();
  if ((nom_elem=="Tetra_VEF")||(nom_elem=="Tri_VEF"))
    istetra=1;
 
  const DoubleVect& porosite_elem = equation().milieu().porosite_elem();
  int ncomp_ch_transporte=(transporte_face.nb_dim() == 1?1:transporte_face.dimension(1));
 
  // Special treatment for periodic faces
  DoubleTrav resu_perio;
  int nb_faces_perio = 0;
  for (int n_bord=0; n_bord<nb_bord; n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
      if (sub_type(Periodique,la_cl.valeur()))
        {
          const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
          nb_faces_perio+=le_bord.nb_faces();
        }
    }
  if (nb_faces_perio>0)
    {
      resu_perio.resize(nb_faces_perio,ncomp_ch_transporte);
      nb_faces_perio=0;
      for (int n_bord=0; n_bord<nb_bord; n_bord++)
        {
          const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
          if (sub_type(Periodique,la_cl.valeur()))
            {
              const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
              int num1 = le_bord.num_premiere_face();
              int num2 = num1 + le_bord.nb_faces();
              CDoubleTabView resu_v = resu.view_ro();
              DoubleTabView resu_perio_v = resu_perio.view_wo();
              Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(num1, num2), KOKKOS_LAMBDA(const int num_face)
              {
                int ind_face = nb_faces_perio + num_face - num1;
                for (int comp = 0; comp < ncomp_ch_transporte; comp++)
                  resu_perio_v(ind_face, comp) = resu_v(num_face, comp);
              });
              end_gpu_timer(__KERNEL_NAME__);
              nb_faces_perio+=num2-num1;
            }
        }
    }
 
  DoubleTab tab_gradient; // Can point to gradient_elem or gradient_face depending on the scheme
  if(type_op==centre || type_op==muscl)
    {
      // Gradient array based on gradient_elem depending on the scheme
      if (gradient_elem_.size_array() == 0) gradient_elem_.resize(nb_elem_tot, ncomp_ch_transporte, dimension);  // (du/dx du/dy dv/dx dv/dy) per element
      Champ_P1NC::calcul_gradient(transporte_face, gradient_elem_, domaine_Cl_VEF);
 
      if (type_op == centre)
        {
          tab_gradient.ref(gradient_elem_);
        }
      else if (type_op == muscl)
        {
          int cas = 1;
          if (type_lim_int == type_lim_minmod) cas = 1;
          if (type_lim_int == type_lim_vanleer) cas = 2;
          if (type_lim_int == type_lim_vanalbada) cas = 3;
          if (type_lim_int == type_lim_chakravarthy) cas = 4;
          if (type_lim_int == type_lim_superbee) cas = 5;
          //  apply the limiter
          if (!gradient_face_.get_md_vector())
            {
              gradient_face_.resize(0, ncomp_ch_transporte, dimension);     // (du/dx du/dy dv/dx dv/dy) per face
              domaine_VEF.creer_tableau_faces(gradient_face_);
            }
          for (int n_bord = 0; n_bord < nb_bord; n_bord++)
            {
              const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
              const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
              int num1 = le_bord.num_premiere_face();
              int num2 = num1 + le_bord.nb_faces();
              if (sub_type(Symetrie, la_cl.valeur()))
                {
                  int ordre = ordre_;
                  int dim = dimension;
                  CIntTabView face_voisins = domaine_VEF.face_voisins().view_ro();
                  CDoubleTabView face_normale = domaine_VEF.face_normales().view_ro();
                  CDoubleTabView3 gradient_elem = gradient_elem_.view_ro<3>();
                  DoubleTabView3 gradient_face = gradient_face_.view_wo<3>();
                  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(num1, num2), KOKKOS_LAMBDA(const int fac)
                  {
                    int elem1 = face_voisins(fac, 0);
                    for (int comp0 = 0; comp0 < ncomp_ch_transporte; comp0++)
                      for (int i = 0; i < dim; i++)
                        gradient_face(fac, comp0, i) = gradient_elem(elem1, comp0, i);
 
                    if (ordre == 3)
                      {
                        // Remove the normal component (could be done for other schemes too...)
                        // but for now, we do not want to change the default behavior of muscl...
                        for (int comp0 = 0; comp0 < ncomp_ch_transporte; comp0++)
                          for (int i = 0; i < dim; i++)
                            {
                              double carre_surface = 0;
                              double tmp = 0;
                              for (int j = 0; j < dim; j++)
                                {
                                  double ndS = face_normale(fac, j);
                                  carre_surface += ndS * ndS;
                                  tmp += gradient_face(fac, comp0, j) * ndS;
                                }
                              gradient_face(fac, comp0, i) -= tmp * face_normale(fac, i) / carre_surface;
                            }
                      }
                  });
                  end_gpu_timer(__KERNEL_NAME__);
                }
            }
          tab_gradient.ref(gradient_face_);
          int dim = dimension;
          int ordre = ordre_;
          CIntTabView face_voisins = domaine_VEF.face_voisins().view_ro();
          CIntArrView traitement_pres_bord = traitement_pres_bord_.view_ro();
          CIntArrView faces_doubles = domaine_VEF.faces_doubles().view_ro();
          CDoubleTabView3 gradient_elem = gradient_elem_.view_ro<3>();
          DoubleTabView3 gradient = tab_gradient.view_wo<3>();
#ifdef TRUST_USE_GPU
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_2D({0,0}, {nb_faces, ncomp_ch_transporte}), KOKKOS_LAMBDA(const int fac, const int comp0)
#else
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(0, nb_faces), KOKKOS_LAMBDA(const int fac)
#endif
          {
            if (faces_doubles[fac] /* face perio */ || fac>=premiere_face_int /* face interne */)
              {
                int elem1 = face_voisins(fac, 0);
                int elem2 = face_voisins(fac, 1);
                int limiteur = cas;
                if (ordre == 3 && (traitement_pres_bord[elem1] || traitement_pres_bord[elem2]))
                  limiteur = 1;
#ifndef TRUST_USE_GPU
                for (int comp0 = 0; comp0 < ncomp_ch_transporte; comp0++)
#endif
                  for (int i = 0; i < dim; i++)
                    {
                      double grad1 = gradient_elem(elem1, comp0, i);
                      double grad2 = gradient_elem(elem2, comp0, i);
                      gradient(fac, comp0, i) = FCT_LIMITEUR(grad1, grad2, limiteur);
                    }
              }
          });
          end_gpu_timer(__KERNEL_NAME__);
          if(nom_elem=="Tetra_VEF")
            {
              /* no echange espace virtuel here, it will be done later in a async way for tetra kernel*/
            }
          else
            {
              tab_gradient.echange_espace_virtuel();// Cannot be removed. Keeping the Kernel on the CPU does not help.
            }
        }// end if(type_op==muscl)
    }
 
  // Dimensioning the array of convective fluxes at the boundary of the computational domain
  DoubleTab& flux_b = flux_bords_;
  int nb_faces_bord=domaine_VEF.nb_faces_bord();
  flux_b.resize(nb_faces_bord,ncomp_ch_transporte);
  // No, lock with DomainFlow test case:
  //if (flux_b.dimension(0)!=nb_faces_bord) flux_b.resize(nb_faces_bord,ncomp_ch_transporte);
  flux_b = 0.;
 
  const IntTab& KEL=type_elemvef.KEL();
  const DoubleTab& xv=domaine_VEF.xv();
  const DoubleTab& coord_sommets=domaine.coord_sommets();
 
 
  // Loop or not depending on the value of alpha (only at order 3 for now)
  // If alpha=1, the loop reduces to a single pass with the chosen scheme (muscl, amont, centre)
  // If alpha<1, the loop consists of 2 passes:
  //                         -the first with the chosen scheme and a weighting of alpha
  //                         -the second with the centre scheme and a weighting of 1-alpha
  double alpha = alpha_;
  int nombre_passes = (alpha==1 ? 1 : 2);
  for (int passe=1; passe<=nombre_passes; passe++)
    {
      type_operateur type_op_boucle = type_op;
      if (passe==2)
        {
          type_op_boucle = centre;
          tab_gradient.ref(gradient_elem_);
        }
      // Non-standard polyhedra are sorted in 2 groups in Domaine_VEF:
      //  - boundary and joint polyhedra
      //  - boundary and non-joint polyhedra
      // Polyhedra are processed in the order they appear
      // in the domain
      // loop over polyhedra
      if(nom_elem=="Tetra_VEF")
        {
 
          FluxTetraKernelData kernel_data;
 
          // Scalars
          kernel_data.nb_elem = nb_elem;
          kernel_data.nb_elem_tot = nb_elem_tot;
          kernel_data.nb_faces = nb_faces;
          kernel_data.nb_faces_bord = nb_faces_bord;
          kernel_data.nfa7 = nfa7;
          kernel_data.ncomp_ch_transporte = ncomp_ch_transporte;
          kernel_data.alpha = alpha;
          kernel_data.marq = marq;
          kernel_data.option_calcul_flux_en_un_point = option_calcul_flux_en_un_point;
          kernel_data.option_appliquer_cl_dirichlet = option_appliquer_cl_dirichlet;
          kernel_data.isMuscl = (type_op_boucle == muscl);
          kernel_data.isAmont = (type_op_boucle == amont);
          // Views
          kernel_data.rang_elem_non_std_v = rang_elem_non_std.view_ro();
          kernel_data.elem_faces_v = elem_faces.view_ro();
          kernel_data.porosite_face_v = porosite_face.view_ro();
          kernel_data.porosite_elem_v = porosite_elem.view_ro();
          kernel_data.coord_sommets_v = coord_sommets.view_ro();
          kernel_data.les_elems_v = les_elems.view_ro();
          kernel_data.facette_normales_v = facette_normales.view_ro<3>();
          kernel_data.est_une_face_de_dirichlet_v = est_une_face_de_dirichlet_.view_ro();
          kernel_data.xp_v = domaine_VEF.xp().view_ro();
          kernel_data.xv_v = xv.view_ro();
          kernel_data.type_elem_Cl_v = type_elem_Cl_.view_ro();
          kernel_data.traitement_pres_bord_v = traitement_pres_bord_.view_ro();
          kernel_data.KEL_v = type_elemvef.KEL().view_ro();
          kernel_data.normales_facettes_Cl_v = normales_facettes_Cl.view_ro<3>();
          kernel_data.vecteur_face_facette_Cl_v = vecteur_face_facette_Cl.view_ro<4>();
          kernel_data.vitesse_v = la_vitesse.valeurs().view_ro();
          kernel_data.transporte_face_v = transporte_face.view_ro();
          kernel_data.gradient_v = tab_gradient.view_ro<3>();
          // Outputs
          kernel_data.resu_v = resu.view_rw();
          kernel_data.flux_b_v = flux_b.view_rw();
 
          if (type_op_boucle == muscl)
            {
              if (ordre_==1)
                {
                  //Wrapper to hide the async complexity
                  compute_flux_tetra_kernel</* ordre */ 1, /*muscl*/ true>(kernel_data, tab_gradient);
                }
 
              else if (ordre_==2)
                {
                  compute_flux_tetra_kernel</* ordre */ 2, /*muscl*/ true>(kernel_data,  tab_gradient);
                }
              else if (ordre_==3)
                {
                  compute_flux_tetra_kernel</* ordre */ 3, /*muscl*/ true>(kernel_data, tab_gradient);
                }
            }
          else
            {
              if (ordre_==1)
                {
                  compute_flux_tetra_kernel</* ordre */ 1, /*muscl*/ false>(kernel_data, tab_gradient);
                }
              else if (ordre_==2)
                {
                  compute_flux_tetra_kernel</* ordre */ 2, /*muscl*/ false>(kernel_data, tab_gradient);
                }
              else if (ordre_==3)
                {
                  compute_flux_tetra_kernel</* ordre */ 3, /*muscl*/ false>(kernel_data, tab_gradient);
                }
            }
        }
      else
        {
          // Non tetra (tri, quad, hexa)
          const DoubleTab& vecteur_face_facette = ref_cast_non_const(Domaine_VEF,domaine_VEF).vecteur_face_facette();
          ArrOfInt face(nfac);
          ArrOfDouble vs(dimension);
          ArrOfDouble vc(dimension);
          ArrOfDouble cc(dimension);
          DoubleVect xc(dimension);
          DoubleTab vsom(nsom,dimension);
          DoubleTab xsom(nsom,dimension);
          for (int poly=0; poly<nb_elem_tot; poly++)
            {
              int contrib = 0;
              // compute the face indices of the polyhedron
              for (int face_adj=0; face_adj<nfac; face_adj++)
                {
                  int face_ = elem_faces(poly,face_adj);
                  face(face_adj)= face_;
                  if (face_<nb_faces) contrib=1; // A real face on the virtual element
                }
              //
              if (contrib)
                {
                  int calcul_flux_en_un_point = (ordre_ != 3) && (ordre_==1 || traitement_pres_bord_(poly));
                  for (int j=0; j<dimension; j++)
                    {
                      vs(j) = velocity_tab(face(0),j)*porosite_face(face(0));
                      for (int i=1; i<nfac; i++)
                        vs(j)+= velocity_tab(face(i),j)*porosite_face(face(i));
                    }
                  // compute velocity at the vertices of the polyhedra
                  // Shape functions implemented in Champs_P1_impl or Champs_Q1_impl will be used
                  if (istetra==1)
                    {
                      for (int i=0; i<nsom; i++)
                        for (int j=0; j<dimension; j++)
                          vsom(i,j) = (vs(j) - dimension*velocity_tab(face(i),j)*porosite_face(face(i)));
                    }
                  else
                    {
                      // to be valid with hexahedra (slower to compute...)
                      for (int j=0; j<nsom; j++)
                        {
                          int num_som = les_elems(poly,j);
                          for (int ncomp=0; ncomp<dimension; ncomp++)
                            vsom(j,ncomp) = la_vitesse.valeur_a_sommet_compo(num_som,poly,ncomp);
                        }
                    }
 
                  // Determine the BC type based on the rank
                  int rang = rang_elem_non_std(poly);
                  int itypcl = (rang==-1 ? 0 : domaine_Cl_VEF.type_elem_Cl(rang));
 
                  // compute vc (at the intersection of the 3 facets); vc, vs, vsom proportional to porosity
                  type_elemvef.calcul_vc(face, vc, vs, vsom, la_vitesse, itypcl, porosite_face);
 
                  // compute xc (at the intersection of the 3 facets), needed for muscl3
                  if (ordre_==3)
                    {
                      int idirichlet;
                      int n1,n2,n3;
                      for (int i=0; i<nsom; i++)
                        for (int j=0; j<dimension; j++)
                          xsom(i,j) = coord_sommets(les_elems(poly,i),j);
                      type_elemvef.calcul_xg(xc,xsom,itypcl,idirichlet,n1,n2,n3);
                    }
 
                  // Porosity handling
                  if (marq==0)
                    {
                      double coeff=1./porosite_elem(poly);
                      vsom*=coeff;
                      vc*=coeff;
                    }
                  // Loop over the facets of the non-standard polyhedron:
                  for (int fa7=0; fa7<nfa7; fa7++)
                    {
                      int num10 = face(KEL(0,fa7));
                      int num20 = face(KEL(1,fa7));
                      // facet normals
                      if (rang==-1)
                        for (int i=0; i<dimension; i++)
                          cc[i] = facette_normales(poly, fa7, i);
                      else
                        for (int i=0; i<dimension; i++)
                          cc[i] = normales_facettes_Cl(rang,fa7,i);
 
                      // Compute velocities at C,S,S2 (the 3 endpoints of fa7) and M (the center of fa7)
                      double psc_c=0,psc_s=0,psc_m,psc_s2=0;
                      if (dimension==2)
                        {
                          for (int i=0; i<2; i++)
                            {
                              psc_c+=vc[i]*cc[i];
                              psc_s+=vsom(KEL(2,fa7),i)*cc[i];
                            }
                          psc_m=(psc_c+psc_s)/2.;
                        }
                      else
                        {
                          for (int i=0; i<3; i++)
                            {
                              psc_c+=vc[i]*cc[i];
                              psc_s+=vsom(KEL(2,fa7),i)*cc[i];
                              psc_s2+=vsom(KEL(3,fa7),i)*cc[i];
                            }
                          psc_m=(psc_c+psc_s+psc_s2)/3.;
                        }
                      // Apply Dirichlet BCs if num1 or num2 is a face with a Dirichlet BC:
                      // in that case fa7 coincides with face num1 or num2 -> C is at the center of the face
                      int appliquer_cl_dirichlet=0;
                      if (option_appliquer_cl_dirichlet)
                        if (est_une_face_de_dirichlet_(num10) || est_une_face_de_dirichlet_(num20))
                          {
                            appliquer_cl_dirichlet = 1;
                            psc_m = psc_c;
                          }
 
                      // Determine the upwind face for point M
                      int face_amont_m,dir;
                      if (psc_m >= 0)
                        {
                          face_amont_m = num10;
                          dir=0;
                        }
                      else
                        {
                          face_amont_m = num20;
                          dir=1;
                        }
                      // Determine the upwind faces for points C,S,S2
                      int face_amont_c=face_amont_m;
                      int face_amont_s=face_amont_m;
                      int face_amont_s2=face_amont_m;
                      if (type_op_boucle==muscl && ordre_==3)
                        {
                          face_amont_c  = (psc_c >= 0)  ? num10 : num20;
                          face_amont_s  = (psc_s >= 0)  ? num10 : num20;
                          face_amont_s2 = (psc_s2 >= 0) ? num10 : num20;
                        }
                      // gradient at element items (centre scheme) or at face items (muscl schemes)
                      int item_m=poly;
                      int item_c=poly;
                      int item_s=poly;
                      int item_s2=poly;
                      if (type_op_boucle==muscl)
                        {
                          item_m = face_amont_m;
                          item_c = face_amont_c;
                          item_s = face_amont_s;
                          item_s2 = face_amont_s2;
                        }
 
                      for (int comp0=0; comp0<ncomp_ch_transporte; comp0++)
                        {
                          double flux;
                          double inco_m = (ncomp_ch_transporte==1?transporte_face(face_amont_m):transporte_face(face_amont_m,comp0));
                          if (type_op_boucle==amont || appliquer_cl_dirichlet)
                            {
                              flux = inco_m*psc_m;
                            }
                          else // muscl or centre
                            {
                              // Compute the unknown at the centre M of fa7
                              if (rang==-1)
                                for (int j=0; j<dimension; j++)
                                  inco_m+= tab_gradient(item_m,comp0,j)*vecteur_face_facette(poly,fa7,j,dir);
                              else
                                for (int j=0; j<dimension; j++)
                                  inco_m+= tab_gradient(item_m,comp0,j)*vecteur_face_facette_Cl(rang,fa7,j,dir);
 
                              // Compute the unknown at vertex S, a first endpoint of fa7
                              double inco_s = (ncomp_ch_transporte==1?transporte_face(face_amont_s):transporte_face(face_amont_s,comp0));
                              int sommet_s = les_elems(poly,KEL(2,fa7));
                              for (int j=0; j<dimension; j++)
                                inco_s+= tab_gradient(item_s,comp0,j)*(-xv(face_amont_s,j)+coord_sommets(sommet_s,j));
 
                              // Compute the unknown at vertex S2, the last endpoint of fa7 in 3D
                              double inco_s2=0;
                              if (dimension==3)
                                {
                                  inco_s2 = (ncomp_ch_transporte==1?transporte_face(face_amont_s2):transporte_face(face_amont_s2,comp0));
                                  int sommet_s2 = les_elems(poly,KEL(3,fa7));
                                  for (int j=0; j<dimension; j++)
                                    inco_s2+= tab_gradient(item_s2,comp0,j)*(-xv(face_amont_s2,j)+coord_sommets(sommet_s2,j));
                                }
 
                              // Compute the unknown at C, another endpoint of fa7, intersection with the other fa7
                              // of the polyhedron. C=G centre of the polyhedron if the volume is not extended
                              // xc given by elemvef.calcul_xg()
                              double inco_c;
                              if (ordre_==3)
                                {
                                  inco_c = (ncomp_ch_transporte==1?transporte_face(face_amont_c):transporte_face(face_amont_c,comp0));
                                  for (int j=0; j<dimension; j++)
                                    inco_c+= tab_gradient(item_c,comp0,j)*(-xv(face_amont_c,j)+xc(j));
                                }
                              else
                                {
                                  inco_c = dimension*inco_m-inco_s-inco_s2;
                                }
 
                              // Compute the flux at 1 point
                              if (calcul_flux_en_un_point || option_calcul_flux_en_un_point)
                                {
                                  flux = inco_m*psc_m;
                                }
                              else
                                {
                                  // Compute the flux at 3 points
                                  flux = (dimension==2) ? (inco_c*psc_c + inco_s*psc_s + 4*inco_m*psc_m)/6
                                         : (inco_c*psc_c + inco_s*psc_s + inco_s2*psc_s2 + 9*inco_m*psc_m)/12;
                                }
                            }
 
                          // Weighting by coefficient alpha
                          flux*=alpha;
 
                          if (ncomp_ch_transporte == 1)
                            {
                              resu(num10) -= flux;
                              resu(num20) += flux;
                              if (num10<nb_faces_bord) flux_b(num10,0) += flux;
                              if (num20<nb_faces_bord) flux_b(num20,0) -= flux;
                            }
                          else
                            {
                              resu(num10,comp0) -= flux;
                              resu(num20,comp0) += flux;
                              if (num10<nb_faces_bord) flux_b(num10,comp0) += flux;
                              if (num20<nb_faces_bord) flux_b(num20,comp0) -= flux;
                            }
 
                        }// loop over components
                    } // end of loop over facets
                }
            } // end of loop
        }
      alpha = 1 - alpha;
    } // end of loop
 
  nb_faces_perio = 0;
  // Loop over the boundaries to process the boundary conditions
  // a convection term is taken into account for
  // Neumann_sortie_libre boundary conditions only
  for (int n_bord=0; n_bord<nb_bord; n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
 
      if (sub_type(Neumann_sortie_libre,la_cl.valeur()))
        {
          const Neumann_sortie_libre& la_sortie_libre = ref_cast(Neumann_sortie_libre, la_cl.valeur());
          const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
          int num1 = le_bord.num_premiere_face();
          int num2 = num1 + le_bord.nb_faces();
          int dim = Objet_U::dimension;
          CDoubleTabView face_normale = domaine_VEF.face_normales().view_ro();
          CDoubleTabView val_ext = la_sortie_libre.val_ext().view_ro();
          CDoubleTabView transporte_face_v = transporte_face.view_ro();
          CDoubleTabView vitesse_face_v = vitesse_face.view_ro();
          DoubleTabView flux_b_v = flux_b.view_wo();
          DoubleTabView resu_v = resu.view_rw();
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(num1, num2), KOKKOS_LAMBDA(const int num_face)
          {
            double psc =0;
            for (int i=0; i<dim; i++)
              psc += vitesse_face_v(num_face,i)*face_normale(num_face,i);
            for (int i=0; i<ncomp_ch_transporte; i++)
              {
                double val = psc > 0 ? transporte_face_v(num_face,i) : val_ext(num_face-num1,i);
                resu_v(num_face,i) -= psc*val;
                flux_b_v(num_face,i) = -psc*val;
              }
          });
          end_gpu_timer(__KERNEL_NAME__);
        }
      else if (sub_type(Periodique,la_cl.valeur()))
        {
          const Periodique& la_cl_perio = ref_cast(Periodique, la_cl.valeur());
          const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
          int num1 = le_bord.num_premiere_face();
          int num2 = num1 + le_bord.nb_faces();
          CIntArrView face_associee = la_cl_perio.face_associee().view_ro();
          CDoubleTabView resu_perio_v = resu_perio.view_ro();
          DoubleTabView resu_v = resu.view_rw();
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(num1, num1+le_bord.nb_faces()/2), KOKKOS_LAMBDA(const int num_face)
          {
            int ind_face = num_face - num1;
            int voisine = face_associee(ind_face) + num1;
            for (int comp=0; comp<ncomp_ch_transporte; comp++)
              {
                double diff1 = resu_v(num_face,comp) - resu_perio_v(nb_faces_perio+ind_face,comp);
                double diff2 = resu_v(voisine,comp)  - resu_perio_v(nb_faces_perio+ind_face+voisine-num_face,comp);
                resu_v(voisine,comp)  += diff1;
                resu_v(num_face,comp) += diff2;
                /* Should not add to flux_b here, it is already computed above
                   flux_b(voisine,comp) += diff1;
                   flux_b(num_face,comp) += diff2; */
              }
          });
          end_gpu_timer(__KERNEL_NAME__);
          nb_faces_perio+=num2-num1;
        }
    }
  modifier_flux(*this);
  return resu;
}
 
// methods retrieved from the old OpVEFFaAmont
KOKKOS_INLINE_FUNCTION void convbisimplicite_dec(const double psc_, const int num1, const int num2,
                                                 CDoubleTabView transporte, const int ncomp,
                                                 Matrice_Morse_View matrice,
                                                 CDoubleArrView porosite_face, int phi_u_transportant)
{
  double psc=psc_;
  if (!phi_u_transportant)
    psc*=porosite_face(psc_>=0 ? num1 : num2);
  for (int comp=0; comp<ncomp; comp++)
    {
      int n1=num1*ncomp+comp;
      int n2=num2*ncomp+comp;
      int n = psc_>=0 ? n1 : n2;
      matrice.atomic_add(n1,n,psc);
      matrice.atomic_add(n2,n,-psc);
    }
}
 
void Op_Conv_VEF_Face::ajouter_contribution(const DoubleTab& tab_transporte, Matrice_Morse& matrice_morse) const
{
  const Champ_Inc_base& v = vitesse();
  ajouter_contribution_gen(tab_transporte, v, matrice_morse);
}
 
void Op_Conv_VEF_Face::ajouter_contribution_gen(const DoubleTab& tab_transporte, const Champ_Inc_base& la_vitesse, Matrice_Morse& matrice_morse) const
{
  modifier_matrice_pour_periodique_avant_contribuer(matrice_morse,equation());
  const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
  const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
 
  const Domaine& domaine = domaine_VEF.domaine();
  const Elem_VEF_base& type_elem = domaine_VEF.type_elem();
  const int nfa7 = type_elem.nb_facette();
  const int nb_elem_tot = domaine_VEF.nb_elem_tot();
 
  int nfac = domaine.nb_faces_elem();
  int nsom = domaine.nb_som_elem();
 
  // For convection treatment, we distinguish standard polyhedra
  // which do not "see" boundary conditions, from non-standard
  // polyhedra which have at least one boundary face.
  // A standard polyhedron has n facets on which the convection
  // scheme is applied.
  // For a non-standard polyhedron with Dirichlet boundary conditions,
  // some facets are carried by boundary faces.
  // In short, for a polyhedron, convection treatment depends
  // on the type (triangle, tetrahedron ...) and the number of Dirichlet faces.
 
  int ncomp_ch_transporte = tab_transporte.nb_dim() == 1 ? 1 : tab_transporte.dimension(1);
 
  int marq=phi_u_transportant(equation());
 
  DoubleTab vitesse_face_;
  // either transporte=phi*transporte_ and vitesse_face=vitesse_
  // or transporte=transporte_ and vitesse_face=phi*vitesse_
  // depending on whether transport is done with phi*u or with u.
  const DoubleTab& tab_velocity_tab = la_vitesse.valeurs();
  const DoubleVect& tab_porosite_face = equation().milieu().porosite_face();
  const DoubleTab& tab_vitesse_face=modif_par_porosite_si_flag(tab_velocity_tab,vitesse_face_,marq,tab_porosite_face);
  const Elem_VEF_base& type_elemvef = domaine_VEF.type_elem();
  int istetra=0;
  Nom nom_elem=type_elemvef.que_suis_je();
  if ((nom_elem=="Tetra_VEF")||(nom_elem=="Tri_VEF"))
    istetra=1;
 
  // Non-standard polyhedra are grouped in 2 categories in Domaine_VEF:
  //  - boundary and joint polyhedra
  //  - boundary and non-joint polyhedra
  // Process polyhedra following the order in which they appear
  // in the domain
 
  // loop over polyhedra
  int phi_u_transportant_yes=phi_u_transportant(equation());
  int dim = Objet_U::dimension;
  CIntTabView KEL = type_elemvef.KEL().view_ro();
  CIntArrView rang_elem_non_std = domaine_VEF.rang_elem_non_std().view_ro();
  CIntArrView type_elem_Cl = domaine_Cl_VEF.type_elem_Cl().view_ro();
  CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
  CDoubleTabView3 normales_facettes_Cl = domaine_Cl_VEF.normales_facettes_Cl().view_ro<3>();
  CDoubleTabView3 facette_normales = domaine_VEF.facette_normales().view_ro<3>();
  CDoubleArrView porosite_face = tab_porosite_face.view_ro();
  CDoubleArrView porosite_elem = equation().milieu().porosite_elem().view_ro();
  CDoubleTabView vitesse = la_vitesse.valeurs().view_ro();
  CDoubleTabView velocity_tab = tab_velocity_tab.view_ro();
  CDoubleTabView transporte = tab_transporte.view_ro();
  Matrice_Morse_View matrice;
  matrice.set(matrice_morse);
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(0, nb_elem_tot), KOKKOS_LAMBDA(const int poly)
  {
    int rang = rang_elem_non_std(poly);
    int itypcl= rang==-1 ? 0 : type_elem_Cl(rang);
 
    // compute the face indices of the polyhedron
    int face[4];
    for (int face_adj=0; face_adj<nfac; face_adj++)
      face[face_adj] = elem_faces(poly,face_adj);
 
    double vs[3];
    for (int j=0; j<dim; j++)
      {
        vs[j] = velocity_tab(face[0],j)*porosite_face[face[0]];
        for (int i=1; i<nfac; i++)
          vs[j] += velocity_tab(face[i],j)*porosite_face[face[i]];
      }
    // compute velocity at the vertices of the polyhedra
    // Shape functions implemented in Champs_P1_impl or Champs_Q1_impl will be used
    double vsom[12];
    if (istetra==1)
      {
        for (int i=0; i<nsom; i++)
          for (int j=0; j<dim; j++)
            vsom[i * dim + j] = (vs[j] - dim*velocity_tab(face[i],j)*porosite_face[face[i]]);
      }
    else
      {
        // to make this valid for hexa elements (it is slower to compute...)
        Process::Kokkos_exit("Implicit scheme for VEF hexa ? Not coded. Even not tested.");
        /*
        for (int i=0; i<nsom; i++)
          {
            int num_som = les_elems(poly,i);
            for (int j=0; j<dim; j++)
              vsom(i,j) = la_vitesse.valeur_a_sommet_compo(num_som,poly,j);
          } */
      }
 
    // compute vc (at the intersection of the 3 facets); vc, vs, vsom proportional to porosity
    double vc[3];
    if (dim==3)
      calcul_vc_tetra_views(face,vc,vs,vsom,vitesse,itypcl,porosite_face);
    else
      calcul_vc_tri_views(face,vc,vs,vsom,vitesse,itypcl,porosite_face);
    if (marq==0)
      {
        const double porosite_poly=porosite_elem(poly);
        for (int j=0; j<dim; j++)
          {
            vs[j]/=porosite_poly;
            vc[j]/=porosite_poly;
            for (int i=0; i<nsom; i++)
              vsom[i * dim + j]/=porosite_poly;
          }
      }
    // Loop over the facets of the non-standard polyhedron:
    double cc[3] = { 0., 0., 0. };
    for (int fa7=0; fa7<nfa7; fa7++)
      {
        // normals to the facets
        for (int i=0; i<dim; i++)
          cc[i] = rang==-1 ? facette_normales(poly,fa7,i) : normales_facettes_Cl(rang,fa7,i);
 
        // Apply the convection scheme to each vertex of the facet
        double psc_c=0,psc_s=0,psc_s2=0;
        for (int i=0; i<dim; i++)
          {
            psc_c += vc[i]*cc[i];
            psc_s += vsom[KEL(2,fa7) * dim + i]*cc[i];
            psc_s2 += dim==3 ? vsom[KEL(3,fa7) * dim + i]*cc[i] : 0;
          }
        double psc_m=(psc_c+psc_s+psc_s2)/dim;
        int num10 = face[KEL(0,fa7)];
        int num20 = face[KEL(1,fa7)];
        convbisimplicite_dec(psc_m,num10,num20,transporte,ncomp_ch_transporte,matrice,porosite_face,phi_u_transportant_yes);
      } // end of loop over facets.
 
  }); // end of loop over polyhedra.
  end_gpu_timer(__KERNEL_NAME__);
 
  // Loop over boundary conditions
  // a convection term is taken into account for
  // Neumann_sortie_libre boundary conditions only
  int nb_bord = domaine_VEF.nb_front_Cl();
  for (int n_bord=0; n_bord<nb_bord; n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
      if (sub_type(Neumann_sortie_libre,la_cl.valeur()))
        {
          const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
          int num1 = le_bord.num_premiere_face();
          int num2 = num1 + le_bord.nb_faces();
          CDoubleTabView face_normales = domaine_VEF.face_normales().view_ro();
          CDoubleTabView vitesse_face = tab_vitesse_face.view_ro();
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(num1, num2), KOKKOS_LAMBDA(const int num_face)
          {
            double psc=0;
            for (int i=0; i<dim; i++)
              psc += vitesse_face(num_face,i)*face_normales(num_face,i);
            if (!phi_u_transportant_yes)
              psc*=porosite_face(num_face);
            if (psc>0)
              {
                for (int j=0; j<ncomp_ch_transporte; j++)
                  {
                    int n0=num_face*ncomp_ch_transporte+j;
                    matrice.atomic_add(n0,n0,psc);
                  }
              }
          });
          end_gpu_timer(__KERNEL_NAME__);
        }
    }
  modifier_matrice_pour_periodique_apres_contribuer(matrice_morse,equation());
}
 
void Op_Conv_VEF_Face::contribue_au_second_membre(DoubleTab& resu ) const
{
  const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
  const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
  const Champ_Inc_base& la_vitesse=vitesse();
  const DoubleTab& face_normales = domaine_VEF.face_normales();
  //const Domaine& domaine = domaine_VEF.domaine();
  //int nfac = domaine.nb_faces_elem();
 
  double psc;
  int i,n_bord;
  int num_face;
  int ncomp;
  if (resu.nb_dim() == 1)
    ncomp=1;
  else
    ncomp= resu.dimension(1);
 
  //ArrOfInt face(nfac);
 
 
  // Special treatment for periodic faces
 
  int nb_faces_perio = 0;
  int voisine;
  double diff1,diff2;
 
  // Loop to count the number of periodic faces
  for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
      if (sub_type(Periodique,la_cl.valeur()))
        {
          const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
          nb_faces_perio += le_bord.nb_faces();
        }
    }
 
  DoubleTab tab;
  if (ncomp == 1)
    tab.resize(nb_faces_perio);
  else
    tab.resize(nb_faces_perio,ncomp);
 
  // Loop to fill tab
  nb_faces_perio=0;
  for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
      if (sub_type(Periodique,la_cl.valeur()))
        {
          const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
          int num1 = le_bord.num_premiere_face();
          int num2 = num1 + le_bord.nb_faces();
          for (num_face=num1; num_face<num2; num_face++)
            {
              if (ncomp == 1)
                tab(nb_faces_perio) = resu(num_face);
              else
                for (int comp=0; comp<ncomp; comp++)
                  tab(nb_faces_perio,comp) = resu(num_face,comp);
              nb_faces_perio++;
            }
        }
    }
  // Loop over the boundaries to process the boundary conditions
  // a convection term is taken into account for
  // Neumann_sortie_libre boundary conditions only
  nb_faces_perio=0;
  for (n_bord=0; n_bord<domaine_VEF.nb_front_Cl(); n_bord++)
    {
 
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
 
      if (sub_type(Neumann_sortie_libre,la_cl.valeur()))
        {
          const Neumann_sortie_libre& la_sortie_libre = ref_cast(Neumann_sortie_libre, la_cl.valeur());
          const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
          int num1 = le_bord.num_premiere_face();
          int num2 = num1 + le_bord.nb_faces();
          for (num_face=num1; num_face<num2; num_face++)
            {
              psc =0;
              for (i=0; i<dimension; i++)
                psc += la_vitesse.valeurs()(num_face,i)*face_normales(num_face,i);
              if (psc>0)
                if (ncomp == 1)
                  resu(num_face) += 0;
                else
                  for (i=0; i<ncomp; i++)
                    resu(num_face,i) += 0;
              else
                {
                  if (ncomp == 1)
                    resu(num_face) -= psc*la_sortie_libre.val_ext(num_face-num1);
                  else
                    for (i=0; i<ncomp; i++)
                      resu(num_face,i) -= psc*la_sortie_libre.val_ext(num_face-num1,i);
                }
            }
        }
      else if (sub_type(Periodique,la_cl.valeur()))
        {
          const Periodique& la_cl_perio = ref_cast(Periodique, la_cl.valeur());
          const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
          int num1 = le_bord.num_premiere_face();
          int num2 = num1 + le_bord.nb_faces();
          ArrOfInt fait(le_bord.nb_faces());
          fait = 0;
          for (num_face=num1; num_face<num2; num_face++)
            {
              if (fait[num_face-num1] == 0)
                {
                  voisine = la_cl_perio.face_associee(num_face-num1) + num1;
 
                  if (ncomp == 1)
                    {
                      diff1 = resu(num_face)-tab(nb_faces_perio);
                      diff2 = resu(voisine)-tab(nb_faces_perio+voisine-num_face);
                      resu(voisine)  += diff1;
                      resu(num_face) += diff2;
                    }
                  else
                    for (int comp=0; comp<ncomp; comp++)
                      {
                        diff1 = resu(num_face,comp)-tab(nb_faces_perio,comp);
                        diff2 = resu(voisine,comp)-tab(nb_faces_perio+voisine-num_face,comp);
                        resu(voisine,comp)  += diff1;
                        resu(num_face,comp) += diff2;
                      }
 
                  fait[num_face-num1]= 1;
                  fait[voisine-num1] = 1;
                }
              nb_faces_perio++;
            }
        }
    }
}
 
void Op_Conv_VEF_Face::remplir_fluent() const
{
  const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
  const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, le_dom_vef.valeur());
  int marq=phi_u_transportant(equation());
  // force computing a time step without "porosity"
  marq=0;
  const Domaine& domaine = domaine_VEF.domaine();
  int nfac = domaine.nb_faces_elem();
  int nsom = domaine.nb_som_elem();
  const int nfa7 = domaine_VEF.type_elem().nb_facette();
  const int nb_elem_tot = domaine_VEF.nb_elem_tot();
 
  // Define the vertex array: (C MALOD 17/07/2007)
 
  // For convection treatment, we distinguish standard polyhedra
  // which do not "see" boundary conditions, from non-standard
  // polyhedra which have at least one boundary face.
  // A standard polyhedron has n facets on which the convection
  // scheme is applied.
  // For a non-standard polyhedron with Dirichlet boundary conditions,
  // some facets are carried by boundary faces.
  // In short, for a polyhedron, convection treatment depends
  // on the type (triangle, tetrahedron ...) and the number of Dirichlet faces.
 
  const Elem_VEF_base& type_elemvef= domaine_VEF.type_elem();
  int istetra=0;
  Nom nom_elem=type_elemvef.que_suis_je();
  if ((nom_elem=="Tetra_VEF")||(nom_elem=="Tri_VEF"))
    istetra=1;
 
  // Reset to zero the array used for
  // the computation of the stability time step
  fluent_ = 0;
 
  // Non-standard polyhedra are grouped in 2 categories in Domaine_VEF:
  //  - boundary and joint polyhedra
  //  - boundary and non-joint polyhedra
  // Process polyhedra following the order in which they appear
  // in the domain
  if (nom_elem=="Tetra_VEF")
    {
      assert(dimension==3);
      int dim = 3;
      CIntArrView rang_elem_non_std = domaine_VEF.rang_elem_non_std().view_ro();
      CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
      CDoubleTabView3 facette_normales = domaine_VEF.facette_normales().view_ro<3>();
      CDoubleTabView3 normales_facettes_Cl = domaine_Cl_VEF.normales_facettes_Cl().view_ro<3>();
      CIntTabView KEL = type_elemvef.KEL().view_ro();
      CDoubleArrView porosite_face = equation().milieu().porosite_face().view_ro();
      CDoubleArrView porosite_elem = equation().milieu().porosite_elem().view_ro();
      CDoubleTabView vitesse_face = vitesse().valeurs().view_ro();
      CIntArrView type_elem_Cl = type_elem_Cl_.view_ro();
      DoubleArrView fluent = fluent_.view_rw();
      // loop over polyhedra
      auto kernel = KOKKOS_LAMBDA(const int poly)
      {
        // compute the face indices of the polyhedron
        int face[4];
        for (int face_adj = 0; face_adj < nfac; face_adj++)
          face[face_adj] = elem_faces(poly, face_adj);
 
        double vs[3];
        for (int j = 0; j < dim; j++)
          {
            vs[j] = 0;
            for (int i = 0; i < nfac; i++)
              vs[j] += vitesse_face(face[i], j) * porosite_face(face[i]);
          }
 
        // compute velocity at the vertices of the tetrahedra
        double vsom[12];
        for (int i = 0; i < 4; i++)
          for (int j = 0; j < 3; j++)
            vsom[i * 3 + j] = (vs[j] - 3 * vitesse_face(face[i], j) * porosite_face(face[i]));
 
        int itypcl = type_elem_Cl(poly);
        // compute vc (at the intersection of the 3 facets); vc, vs, vsom proportional to porosity
        double vc[3];
        calcul_vc_tetra_views(face, vc, vs, vsom, vitesse_face, itypcl, porosite_face);
        if (marq == 0)
          {
            double inv_porosite_poly = 1.0 / porosite_elem(poly);
            for (int j = 0; j < dim; j++)
              {
                vs[j] *= inv_porosite_poly;
                vc[j] *= inv_porosite_poly;
                for (int i = 0; i < nsom; i++)
                  vsom[i * dim + j] *= inv_porosite_poly;
              }
          }
        int rang = rang_elem_non_std(poly);
        // Loop over the facets of the non-standard polyhedron:
        for (int fa7 = 0; fa7 < nfa7; fa7++)
          {
            int num1 = face[KEL(0, fa7)];
            int num2 = face[KEL(1, fa7)];
            // Apply the convection scheme to each vertex of the facet
            double psc_c = 0, psc_s = 0, psc_s2 = 0;
            for (int i = 0; i < dim; i++)
              {
                double cc = rang == -1 ? facette_normales(poly, fa7, i) : normales_facettes_Cl(rang, fa7, i);
                psc_c  += vc[i] * cc;
                psc_s  += vsom[KEL(2, fa7) * dim + i] * cc;
                psc_s2 += vsom[KEL(3, fa7) * dim + i] * cc;
              }
            double psc_m = (psc_c + psc_s + psc_s2) / dim;
 
            int num = (psc_m >= 0 ? num2 : num1);
            Kokkos::atomic_add(&fluent[num], std::abs(psc_m));
          } // end of loop over facets
      };
      Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_elem_tot, kernel);
      end_gpu_timer(__KERNEL_NAME__);
    }
  else
    {
      const DoubleVect& porosite_elem = equation().milieu().porosite_elem();
      const IntTab& KEL = type_elemvef.KEL();
      const DoubleVect& porosite_face = equation().milieu().porosite_face();
      const IntTab& elem_faces = domaine_VEF.elem_faces();
      const auto& facette_normales = domaine_VEF.facette_normales();
      const IntTab& les_elems = domaine.les_elems();
      const IntVect& rang_elem_non_std = domaine_VEF.rang_elem_non_std();
      const DoubleTab& normales_facettes_Cl = domaine_Cl_VEF.normales_facettes_Cl();
      ArrOfInt face(nfac);
      ArrOfDouble vs(dimension);
      ArrOfDouble vc(dimension);
      DoubleTab vsom(nsom,dimension);
      ArrOfDouble cc(dimension);
      const DoubleTab& vitesse_face = vitesse().valeurs();
      // loop over polyhedra
      for (int poly=0; poly<nb_elem_tot; poly++)
        {
          int rang = rang_elem_non_std(poly);
          // For an element that is not on the boundary (rang==-1),
          // check if one of its vertices is on the boundary (vertex array) (C MALOD 17/07/2007)
          int itypcl;
          if (rang==-1)
            itypcl=0;
          else
            itypcl=domaine_Cl_VEF.type_elem_Cl(rang);
 
          // compute the face indices of the polyhedron
          for (int face_adj=0; face_adj<nfac; face_adj++)
            face[face_adj]= elem_faces(poly,face_adj);
 
          for (int j=0; j<dimension; j++)
            {
              vs[j] = vitesse_face(face[0],j)*porosite_face[face[0]];
              for (int i=1; i<nfac; i++)
                vs[j]+= vitesse_face(face[i],j)*porosite_face[face[i]];
            }
          // compute velocity at the vertices of the polyhedra
          // Shape functions implemented in Champs_P1_impl or Champs_Q1_impl will be used
          if (istetra==1)
            {
              for (int i=0; i<nsom; i++)
                for (int j=0; j<dimension; j++)
                  vsom(i,j) = (vs[j] - dimension*vitesse_face(face[i],j)*porosite_face[face[i]]);
            }
          else
            {
              // to be valid with hexahedra (slower to compute...)
              int ncomp;
              for (int j=0; j<nsom; j++)
                {
                  int num_som = les_elems(poly,j);
                  for (ncomp=0; ncomp<dimension; ncomp++)
                    vsom(j,ncomp) = vitesse().valeur_a_sommet_compo(num_som,poly,ncomp);
                }
            }
 
          // compute vc (at the intersection of the 3 facets); vc, vs, vsom proportional to porosity
          type_elemvef.calcul_vc(face,vc,vs,vsom,vitesse(),itypcl,porosite_face);
          if (marq==0)
            {
              double porosite_poly=porosite_elem(poly);
              for (int i=0; i<nsom; i++)
                for (int j=0; j<dimension; j++)
                  vsom(i,j)/=porosite_poly;
              for (int j=0; j<dimension; j++)
                {
                  vs[j]/= porosite_poly;
                  vc[j]/=porosite_poly;
                }
            }
          // Loop over the facets of the non-standard polyhedron:
          for (int fa7=0; fa7<nfa7; fa7++)
            {
              int num1 = face[KEL(0,fa7)];
              int num2 = face[KEL(1,fa7)];
              // facet normals
              if (rang==-1)
                {
                  for (int i=0; i<dimension; i++)
                    cc[i] = facette_normales(poly, fa7, i);
                }
              else
                {
                  for (int i=0; i<dimension; i++)
                    cc[i] = normales_facettes_Cl(rang,fa7,i);
                }
              // Apply the convection scheme to each vertex of the facet
 
              double psc_c=0,psc_s=0,psc_m,psc_s2=0;
              if (dimension==2)
                {
                  for (int i=0; i<dimension; i++)
                    {
                      psc_c+=vc[i]*cc[i];
                      psc_s+=vsom(KEL(2,fa7),i)*cc[i];
                    }
                  psc_m=(psc_c+psc_s)/2.;
                }
              else
                {
                  for (int i=0; i<dimension; i++)
                    {
                      psc_c+=vc[i]*cc[i];
                      psc_s+=vsom(KEL(2,fa7),i)*cc[i];
                      psc_s2+=vsom(KEL(3,fa7),i)*cc[i];
                    }
                  psc_m=(psc_c+psc_s+psc_s2)/3.;
                }
 
              // int amont,dir;
              if (psc_m >= 0)
                {
                  // amont = num1;
                  fluent_(num2)  += psc_m;
                  //dir=0;
                }
              else
                {
                  //amont = num2;
                  fluent_(num1)  -= psc_m;
                  //dir=1;
                }
 
            } // end of loop over facets
        } // end of loop
    }
 
  // Loop over boundary conditions
  // a convection term is taken into account for
  // Neumann_sortie_libre boundary conditions only
  int nb_bord = domaine_VEF.nb_front_Cl();
  int nb_comp = Objet_U::dimension;
  for (int n_bord=0; n_bord<nb_bord; n_bord++)
    {
      const Cond_lim& la_cl = domaine_Cl_VEF.les_conditions_limites(n_bord);
      if (sub_type(Neumann_sortie_libre,la_cl.valeur()))
        {
          const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
          int num1b = le_bord.num_premiere_face();
          int num2b = num1b + le_bord.nb_faces();
          CDoubleTabView face_normales = domaine_VEF.face_normales().view_ro();
          CDoubleTabView vitesse_face = vitesse().valeurs().view_ro();
          DoubleArrView fluent = fluent_.view_rw();
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                               Kokkos::RangePolicy<>(num1b, num2b), KOKKOS_LAMBDA(
                                 const int num_face)
          {
            double psc = 0;
            for (int i=0; i<nb_comp; i++)
              psc += vitesse_face(num_face,i)*face_normales(num_face,i);
            if (psc>0) { /* Do nothing */}
            else
              fluent(num_face) -= psc;
          });
          end_gpu_timer(__KERNEL_NAME__);
        }
    }
}
 
void Op_Conv_VEF_Face::get_ordre(int& ord) const
{
  ord=ordre_;
}
void Op_Conv_VEF_Face::get_type_lim(Motcle& typelim) const
{
  typelim=type_lim;
}
void Op_Conv_VEF_Face::get_alpha(double& alp) const
{
  alp=alpha_;
}
 
void Op_Conv_VEF_Face::get_type_op(int& typeop) const
{
  typeop=type_op;
}