Op_Div_VEFP1B_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 <Dirichlet_entree_fluide_leaves.h>
#include <Check_espace_virtuel.h>
#include <Op_Div_VEFP1B_Elem.h>
#include <Dirichlet_homogene.h>
#include <Schema_Temps_base.h>
#include <Porosites_champ.h>
#include <Neumann_val_ext.h>
#include <Domaine_Cl_VEF.h>
#include <communications.h>
#include <EcrFicPartage.h>
#include <Probleme_base.h>
#include <VerifierCoin.h>
#include <Domaine_VEF.h>
#include <Periodique.h>
#include <Dirichlet.h>
#include <Symetrie.h>
#include <Domaine.h>
#include <SChaine.h>
#include <Neumann.h>
#include <Debog.h>
#include <Device.h>
 
Implemente_instanciable(Op_Div_VEFP1B_Elem, "Op_Div_VEFPreP1B_P1NC|Op_Div_VEF_P1NC", Operateur_Div_base);
 
Sortie& Op_Div_VEFP1B_Elem::printOn(Sortie& s) const { return s << que_suis_je() ; }
 
Entree& Op_Div_VEFP1B_Elem::readOn(Entree& is) { return is; }
 
static int chercher_arete(int elem, int somi, int somj, const IntTab& elem_aretes, const IntTab& aretes_som)
{
  if (somi > somj)
    {
      int k = somi;
      somi = somj;
      somj = k;
    }
  for (int i_arete = 0; i_arete < 6; i_arete++)
    {
      int arete = elem_aretes(elem, i_arete);
      int som1 = aretes_som(arete, 0);
      if (somi == som1)
        {
          int som2 = aretes_som(arete, 1);
          if (somj == som2)
            return arete;
        }
    }
  return -1;
}
 
static int verifier(const Op_Div_VEFP1B_Elem& op, int& init, const Domaine_VEF& domaine_VEF, const DoubleTab& vit, DoubleTab& div)
{
  init = 1;
  DoubleTab v(vit);
  v = 1;
  DoubleTab r(div);
  r = 0;
  op.ajouter(v, r);
  Cerr << "div(1,..,1) = " << r << finl;
  Debog::verifier("div(1,..,1) =", r);
  const DoubleTab& xv = domaine_VEF.xv();
  v = xv;
  {
    for (int i = 0; i < v.dimension(0); i++)
      v(i, 1) = 0;
  }
  r = 0;
  op.ajouter(v, r);
  Cerr << "div(x,0) = " << r << finl;
  Debog::verifier("div(x,0) =", r);
  return 1;
}
 
void Op_Div_VEFP1B_Elem::associer(const Domaine_dis_base& domaine_dis, const Domaine_Cl_dis_base& domaine_Cl_dis, const Champ_Inc_base&)
{
  le_dom_vef = ref_cast(Domaine_VEF, domaine_dis);
  la_zcl_vef = ref_cast(Domaine_Cl_VEF, domaine_Cl_dis);
}
 
DoubleTab& Op_Div_VEFP1B_Elem::ajouter_elem(const DoubleTab& vit, DoubleTab& div) const
{
  const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, le_dom_vef.valeur());
  assert(domaine_VEF.get_alphaE());
  const Domaine& domaine = domaine_VEF.domaine();
  const DoubleTab& face_normales = domaine_VEF.face_normales();
  const IntTab& elem_faces = domaine_VEF.elem_faces();
  const IntTab& face_voisins = domaine_VEF.face_voisins();
  int nfe = domaine.nb_faces_elem();
  int nb_elem = domaine.nb_elem();
  int dim = Objet_U::dimension;  // Objet_U::dimension can not be read from Kernel.
 
  CIntTabView elem_faces_v = elem_faces.view_ro();
  DoubleTabView div_v = div.view_rw(); // read-write
  if (getenv("TRUST_USE_RANDOM_ACCESS")!=nullptr)
    {
      // Random access
      RandomAccessView<int, 2> face_voisins_v = face_voisins.view_ro();
      RandomAccessView<double, 2> face_normales_v = face_normales.view_ro();
      RandomAccessView<double, 2> vit_v = vit.view_ro();
      Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_elem, KOKKOS_LAMBDA(
                             const int elem)
      {
        double pscf = 0;
        for (int indice = 0; indice < nfe; indice++)
          {
            int face = elem_faces_v(elem, indice);
            int signe = elem == face_voisins_v(face, 0) ? 1 : -1;
            for (int comp = 0; comp < dim; comp++)
              pscf += signe * vit_v(face, comp) * face_normales_v(face, comp);
          }
        div_v(elem, 0) += pscf;
      });
    }
  else
    {
      CIntTabView face_voisins_v = face_voisins.view_ro();
      CDoubleTabView face_normales_v = face_normales.view_ro();
      CDoubleTabView vit_v = vit.view_ro();
      Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_elem, KOKKOS_LAMBDA(
                             const int elem)
      {
        double pscf = 0;
        for (int indice = 0; indice < nfe; indice++)
          {
            int face = elem_faces_v(elem, indice);
            int signe = elem == face_voisins_v(face, 0) ? 1 : -1;
            for (int comp = 0; comp < dim; comp++)
              pscf += signe * vit_v(face, comp) * face_normales_v(face, comp);
          }
        div_v(elem, 0) += pscf;
      });
    }
  end_gpu_timer(__KERNEL_NAME__);
  assert_invalide_items_non_calcules(div);
  return div;
}
 
int find_cl_face(const Domaine& domaine, const int face)
{
  int i;
  const int nb_cl = domaine.nb_front_Cl();
  for (i = 0; i < nb_cl; i++)
    {
      const Frontiere& fr = domaine.frontiere(i);
      const int debut = fr.num_premiere_face();
      if (face >= debut && face < debut + fr.nb_faces())
        return i;
      const ArrOfInt& faces_virt = fr.get_faces_virt();
      const int nb_faces_virt = faces_virt.size_array();
      for (int j = 0; j < nb_faces_virt; j++)
        if (face == faces_virt[j])
          return i;
    }
  return -1;
}
 
KOKKOS_FUNCTION
double calculer_coef_som(int type_elem, int dimension, int& nb_face_diri, int* indice_diri)
{
  if (type_elem == 0)
    {
      nb_face_diri = 0;
      double coeff_som = 1. / (dimension * (dimension + 1));
      return coeff_som;
    }
  else
    {
      if (dimension == 2)
        {
          switch(type_elem)
            {
            case 1:
              nb_face_diri = 1;
              indice_diri[0] = 2;
              break;
            case 2:
              nb_face_diri = 1;
              indice_diri[0] = 1;
              break;
            case 4:
              nb_face_diri = 1;
              indice_diri[0] = 0;
              break;
            case 3:
              nb_face_diri = 2;
              indice_diri[0] = 1;
              indice_diri[1] = 2;
              break;
            case 5:
              nb_face_diri = 2;
              indice_diri[0] = 0;
              indice_diri[1] = 2;
              break;
            case 6:
              nb_face_diri = 2;
              indice_diri[0] = 0;
              indice_diri[1] = 1;
              break;
            default:
              //abort();
              break;
            }
        }
      else if (dimension == 3)
        {
          switch(type_elem)
            {
            case 1:
              nb_face_diri = 1;
              indice_diri[0] = 3;
              break;
            case 2:
              nb_face_diri = 1;
              indice_diri[0] = 2;
              break;
            case 3:
              nb_face_diri = 2;
              indice_diri[0] = 3;
              indice_diri[1] = 2;
              break;
            case 4:
              nb_face_diri = 1;
              indice_diri[0] = 1;
              break;
            case 5:
              nb_face_diri = 2;
              indice_diri[0] = 1;
              indice_diri[1] = 3;
              break;
            case 6:
              nb_face_diri = 2;
              indice_diri[0] = 1;
              indice_diri[1] = 2;
              break;
            case 7:
              nb_face_diri = 3;
              indice_diri[0] = 1;
              indice_diri[1] = 2;
              indice_diri[2] = 3;
              break;
            case 8:
              nb_face_diri = 1;
              indice_diri[0] = 0;
              break;
            case 9:
              nb_face_diri = 2;
              indice_diri[0] = 0;
              indice_diri[1] = 3;
              break;
            case 10:
              nb_face_diri = 2;
              indice_diri[0] = 0;
              indice_diri[1] = 2;
              break;
            case 11:
              nb_face_diri = 3;
              indice_diri[0] = 0;
              indice_diri[1] = 2;
              indice_diri[2] = 3;
              break;
            case 12:
              nb_face_diri = 2;
              indice_diri[0] = 0;
              indice_diri[1] = 1;
              break;
            case 13:
              nb_face_diri = 3;
              indice_diri[0] = 0;
              indice_diri[1] = 1;
              indice_diri[2] = 3;
              break;
            case 14:
              nb_face_diri = 3;
              indice_diri[0] = 0;
              indice_diri[1] = 1;
              indice_diri[2] = 2;
              break;
            default:
#ifndef TRUST_USE_GPU
              abort();
#endif
              break;
            }
        }
      else
        {
#ifndef TRUST_USE_GPU
          abort();
#endif
        }
    }
  double coeff_som = 1. / (dimension * (dimension + 1 - nb_face_diri));
  return coeff_som;
}
 
DoubleTab& Op_Div_VEFP1B_Elem::ajouter_som(const DoubleTab& tab_vit, DoubleTab& tab_div, DoubleTab& tab_flux_b) const
{
  const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, le_dom_vef.valeur());
  assert(domaine_VEF.get_alphaS());
  const Domaine& domaine = domaine_VEF.domaine();
  int nfe = domaine.nb_faces_elem();
  int nb_elem_tot = domaine.nb_elem_tot();
  int nps = domaine_VEF.numero_premier_sommet();
 
  // Initialize constant arrays
  if (!som_initialized_)
    {
      som_initialized_ = true;
      const IntTab& som_elem = domaine.les_elems();
      som_.resize(nb_elem_tot, nfe);
      nb_degres_liberte_.resize(domaine_VEF.domaine().nb_som_tot());
      nb_degres_liberte_ = -1;
      for (int elem = 0; elem < nb_elem_tot; elem++)
        for (int indice = 0; indice < nfe; indice++)
          {
            int som = nps + domaine.get_renum_som_perio(som_elem(elem, indice));
            nb_degres_liberte_(som - nps)++;
            som_(elem, indice) = som;
          }
      corrige_sommets_sans_degre_liberte_ = (mp_min_vect(nb_degres_liberte_) == 0);
    }
 
  int dim = Objet_U::dimension;
  int modif_traitement_diri = domaine_VEF.get_modif_div_face_dirichlet();
  const Domaine_Cl_VEF& zcl = ref_cast(Domaine_Cl_VEF, la_zcl_vef.valeur());
 
  CDoubleTabView face_normales = domaine_VEF.face_normales().view_ro();
  CDoubleTabView vit = tab_vit.view_ro();
  CIntArrView rang_elem_non_std = domaine_VEF.rang_elem_non_std().view_ro();
  CIntArrView type_elem_Cl = zcl.type_elem_Cl().view_ro();
  CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
  CIntTabView face_voisins = domaine_VEF.face_voisins().view_ro();
  CIntTabView som_v = som_.view_ro();
  DoubleArrView div = static_cast<DoubleVect&>(tab_div).view_rw();
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                       range_1D(0, nb_elem_tot),
                       KOKKOS_LAMBDA (const int elem)
  {
    double sigma[3];
    for (int comp = 0; comp < dim; comp++)
      {
        double s = 0;
        for (int indice = 0; indice < nfe; indice++)
          {
            int face = elem_faces(elem, indice);
            s += vit(face, comp);
          }
        sigma[comp] = s;
      }
 
    double coeff_som = 1. / (dim * (dim + 1));
    if (modif_traitement_diri)
      {
        int indice_diri[4];
        int nb_face_diri = 0;
        int rang_elem = (int)rang_elem_non_std(elem);
        int type_elem = rang_elem < 0 ? 0 : (int)type_elem_Cl(rang_elem);
        coeff_som = calculer_coef_som(type_elem, dim, nb_face_diri, indice_diri);
        // subtract the contribution from Dirichlet faces
        for (int fdiri = 0; fdiri < nb_face_diri; fdiri++)
          {
            int indice = indice_diri[fdiri];
            int face = elem_faces(elem,indice);
            for (int comp = 0; comp < dim; comp++)
              sigma[comp] -= vit(face,comp);
          }
      }
 
    for (int indice = 0; indice < nfe; indice++)
      {
        int face = elem_faces(elem,indice);
        double psc = 0;
        for (int comp = 0; comp < dim; comp++)
          psc += sigma[comp] * face_normales(face,comp);
 
        int som = som_v(elem,indice);
        int signe = (elem != face_voisins(face,0) ? -1 : 1);
        Kokkos::atomic_add(&div(som), signe * coeff_som * psc);
      }
  });
  end_gpu_timer(__KERNEL_NAME__);
 
  const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
  const Conds_lim& les_cl = domaine_Cl_VEF.les_conditions_limites();
  int nb_bords = les_cl.size();
  int nb_comp = Objet_U::dimension;
 
  IntArrView nb_degres_liberte = nb_degres_liberte_.view_rw();
  DoubleTabView flux_b = tab_flux_b.view_wo();
 
  for (int n_bord = 0; n_bord < nb_bords; 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());
      CIntArrView num_face = le_bord.num_face().view_ro();
      int nb_faces_bord = le_bord.nb_faces();
      int nb_faces_bord_tot = le_bord.nb_faces_tot();
      assert(le_bord.nb_faces() == domaine_VEF.domaine().frontiere(n_bord).nb_faces());
      if (!sub_type(Periodique, la_cl.valeur()))
        {
          int libre = 1;
          if (sub_type(Dirichlet,la_cl.valeur()) || sub_type(Dirichlet_homogene, la_cl.valeur()) || sub_type(Dirichlet_entree_fluide, la_cl.valeur()) || sub_type(Symetrie, la_cl.valeur()))
            libre = 0;
 
          CIntTabView face_sommets = domaine_VEF.face_sommets().view_ro();
          CIntArrView renum_som_perio = domaine.get_renum_som_perio().view_ro();
 
          // Loop over real and virtual boundary faces
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                               range_1D(0, nb_faces_bord_tot), KOKKOS_LAMBDA(
                                 const int ind_face)
          {
            int face = num_face(ind_face);
            double flux = 0.;
            for (int comp = 0; comp < nb_comp; comp++)
              flux += vit(face, comp) * face_normales(face, comp);
            if (ind_face < nb_faces_bord)
              flux_b(face, 0) = flux;
            flux *= 1. / nb_comp;
            for (int indice = 0; indice < (nfe - 1); indice++)
              {
                int som = renum_som_perio(face_sommets(face, indice));
                Kokkos::atomic_add(&div(nps + som), flux);
                if (libre)
                  Kokkos::atomic_add(&nb_degres_liberte(som), 1);
              }
          });
          end_gpu_timer(__KERNEL_NAME__);
        }
      else
        {
          const Periodique& la_cl_perio = ref_cast(Periodique, la_cl.valeur());
          CIntArrView face_associee = la_cl_perio.face_associee().view_ro();
          // Loop over real and virtual boundary faces
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                               range_1D(0, nb_faces_bord_tot), KOKKOS_LAMBDA(
                                 const int ind_face)
          {
            int face = num_face(ind_face);
            int face_perio = num_face(face_associee(ind_face));
            double flux = 0.;
            double flux_perio = 0.;
            for (int comp = 0; comp < nb_comp; comp++)
              {
                flux += vit(face, comp) * face_normales(face, comp);
                flux_perio += vit(face_perio, comp) * face_normales(face_perio, comp);
              }
            if (ind_face < (nb_faces_bord / 2))
              {
                flux_b(face, 0) = -flux;
                flux_b(face_perio, 0) = flux_perio;
              }
          });
          end_gpu_timer(__KERNEL_NAME__);
        }
    }
  return tab_div;
}
 
DoubleTab& Op_Div_VEFP1B_Elem::ajouter_aretes(const DoubleTab& vit, DoubleTab& div) const
{
  const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, le_dom_vef.valeur());
  assert(domaine_VEF.get_alphaA());
  const Domaine& domaine = domaine_VEF.domaine();
  const DoubleTab& face_normales = domaine_VEF.face_normales();
  const IntTab& som_elem = domaine.les_elems();
  const IntTab& elem_faces = domaine_VEF.elem_faces();
  const IntTab& face_voisins = domaine_VEF.face_voisins();
  const ArrOfInt& renum_arete_perio = domaine_VEF.get_renum_arete_perio();
  const ArrOfInt& ok_arete = domaine_VEF.get_ok_arete();
  int npa = domaine_VEF.numero_premiere_arete();
  int nb_elem_tot = domaine.nb_elem_tot();
  int elem, comp;
 
  const IntTab& aretes_som = domaine_VEF.domaine().aretes_som();
  const IntTab& elem_aretes = domaine_VEF.domaine().elem_aretes();
  for (elem = 0; elem < nb_elem_tot; elem++)
    {
      for (int isom = 0; isom < 3; isom++)
        {
          int somi = som_elem(elem, isom);
          int facei = elem_faces(elem, isom);
          double signei = 1.;
          if (face_voisins(facei, 0) != elem)
            signei = -1.;
          for (int jsom = isom + 1; jsom < 4; jsom++)
            {
              int somj = som_elem(elem, jsom);
              int facej = elem_faces(elem, jsom);
              double signej = 1.;
              if (face_voisins(facej, 0) != elem)
                signej = -1.;
              int arete = renum_arete_perio[chercher_arete(elem, somi, somj, elem_aretes, aretes_som)];
              if (ok_arete[arete])
                {
                  int niinij = 0;
                  for (int ksom = 0; ksom < 4; ksom++)
                    {
                      if ((ksom == isom) || (ksom == jsom))
                        {
                          if (niinij == ksom)
                            niinij++;
                          double psc = 0;
                          for (comp = 0; comp < 3; comp++)
                            psc += (signei * face_normales(facei, comp) + signej * face_normales(facej, comp)) * vit(elem_faces(elem, ksom), comp);
                          div(npa + arete) -= 1. / 15. * psc;
                        }
                      else
                        {
                          if (niinij == ksom)
                            niinij++;
                          while ((niinij == jsom) || (niinij == isom))
                            niinij++;
                          assert(niinij < 4);
                          assert(niinij != isom);
                          assert(niinij != jsom);
                          assert(niinij != ksom);
                          int facek = elem_faces(elem, niinij);
                          double signek = 1.;
                          if (face_voisins(facek, 0) != elem)
                            signek = -1.;
                          double psc = 0;
                          for (comp = 0; comp < 3; comp++)
                            psc += signek * face_normales(facek, comp) * vit(elem_faces(elem, ksom), comp);
                          div(npa + arete) -= 2. / 15. * psc;
                          niinij = ksom;
                        }
                    }
                }
            }
        }
    }
  return div;
}
 
DoubleTab& Op_Div_VEFP1B_Elem::ajouter(const DoubleTab& tab_velocity_tab, DoubleTab& tab_div) const
{
  const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, le_dom_vef.valeur());
  const Domaine_Cl_VEF& domaine_Cl_VEF = la_zcl_vef.valeur();
  // Some checks:
  // The virtual space of velocity_tab must be up to date (test is done if check_enabled==1)
  assert_espace_virtuel_vect(tab_velocity_tab);
  // We don't care about the virtual space of div at input, but we do += on it.
  assert_invalide_items_non_calcules(tab_div, 0.);
 
#ifndef NDEBUG
  // Ensure periodicity is enforced on velocity_tab (see FA814)
  int nb_comp = tab_velocity_tab.dimension(1);
  for (int 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 Periodique& la_cl_perio = ref_cast(Periodique, la_cl.valeur());
          const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
          int nb_faces_bord = le_bord.nb_faces();
          CIntArrView face_associee = la_cl_perio.face_associee().view_ro();
          CIntArrView num_face = le_bord.num_face().view_ro();
          CDoubleTabView velocity_tab = tab_velocity_tab.view_ro();
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(0, nb_faces_bord), KOKKOS_LAMBDA(const int ind_face)
          {
            int ind_face_associee = face_associee(ind_face);
            int face = num_face(ind_face);
            int face_ass = num_face(ind_face_associee);
            // PL: tolerance 1.e-4 because some cases do not satisfy this strictly (slow drift?)
            for (int comp = 0; comp < nb_comp; comp++)
              if (!est_egal(velocity_tab(face, comp), velocity_tab(face_ass, comp), 1.e-4))
                {
#ifndef TRUST_USE_GPU
                  Cerr << "vit1(" << face << "," << comp << ")=" << velocity_tab(face, comp) << finl;
                  Cerr << "vit2(" << face_ass << "," << comp << ")=" << velocity_tab(face_ass, comp) << finl;
                  Cerr << "Delta=" << velocity_tab(face, comp) - velocity_tab(face_ass, comp) << finl;
                  Cerr << "Periodic boundary condition is not correct in Op_Div_VEFP1B_Elem::ajouter" << finl;
                  Cerr << "Contact TRUST support." << finl;
#endif
                  Process::Kokkos_exit("Error in Op_Div_VEFP1B_Elem::ajouter.");
                }
          });  // for face
          end_gpu_timer(__KERNEL_NAME__);
        }  // sub_type Perio
    }
#endif
  const DoubleVect& porosite_face = equation().milieu().porosite_face();
  DoubleTab phi_vitesse_face_;
  const DoubleTab& vit = modif_par_porosite_si_flag(tab_velocity_tab, phi_vitesse_face_, 1, porosite_face);
 
  DoubleTab& flux_b = flux_bords_;
  flux_b.resize(domaine_VEF.nb_faces_bord(), 1);
  flux_b = 0.;
 
  static int init = 1;
  if (!init)
    ::verifier(*this, init, domaine_VEF, vit, tab_div);
  if (domaine_VEF.get_alphaE())
    ajouter_elem(vit, tab_div);
  if (domaine_VEF.get_alphaS())
    ajouter_som(vit, tab_div, flux_b);
  if (domaine_VEF.get_alphaA())
    ajouter_aretes(vit, tab_div);
 
  // correction of div u if vertex pressure is strongly imposed
  if ((domaine_VEF.get_alphaS()) && ((domaine_VEF.get_cl_pression_sommet_faible() == 0)))
    {
      const Conds_lim& les_cl = domaine_Cl_VEF.les_conditions_limites();
      int nps = domaine_VEF.numero_premier_sommet();
      for (const auto &itr : les_cl)
        {
          const Cond_lim& la_cl = itr;
          if (sub_type(Neumann,la_cl.valeur()) || sub_type(Neumann_val_ext, la_cl.valeur()))
            {
              const Front_VF& la_front_dis = ref_cast(Front_VF, la_cl->frontiere_dis());
              int nb_faces = la_front_dis.nb_faces_tot();
              int nsf = 0;
              if (nb_faces != 0)
                nsf = domaine_VEF.face_sommets().dimension(1);
              CIntArrView num_face = la_front_dis.num_face().view_ro();
              CIntTabView faces_sommets = domaine_VEF.face_sommets().view_ro();
              DoubleArrView div = static_cast<DoubleVect&>(tab_div).view_wo();
              Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                                   range_1D(0, nb_faces), KOKKOS_LAMBDA(
                                     const int ind_face)
              {
                int face = num_face(ind_face);
                for (int som = 0; som < nsf; som++)
                  {
                    int som1 = faces_sommets(face, som);
                    div(nps + som1) = 0.;
                  }
              });
              end_gpu_timer(__KERNEL_NAME__);
            }
        }
    }
  if (domaine_VEF.get_alphaS() && corrige_sommets_sans_degre_liberte_)
    degres_liberte();
  //Optimization, not necessary:
  //tab_div.echange_espace_virtuel();
  return tab_div;
}
 
#ifdef VersionP1
double tiers=1./3.;
const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
const IntTab& elem_aretes=domaine_VEF.domaine().elem_aretes();
ArrOfDouble gradK(3);
for(int isom=0; isom<3; isom++)
  {
    int somi = som_elem(elem,isom);
    int facei = elem_faces(elem,isom);
    double signei=1.;
    if (face_voisins(facei,0) != elem)
      signei=-1.;
    for(int jsom=isom+1; jsom<4; jsom++)
      {
        int somj = som_elem(elem,jsom);
        int facej = elem_faces(elem,jsom);
        double signej=1.;
        if (face_voisins(facej,0) != elem)
          signej=-1.;
        int arete = chercher_arete(elem, somi, somj,
                                   elem_aretes, aretes_som);
        for(comp=0; comp<dimension; comp++)
          gradK(comp) = (signei*face_normales(facei,comp)+
                         signej*face_normales(facej,comp));
        psc=0;
        for(comp=0; comp<dimension; comp++)
          psc+=tiers*(
                 -signei*face_normales(facei,comp)*vit(facei,comp)
                 -signej*face_normales(facej,comp)*vit(facej,comp)
                 +.5*gradK(comp)*sigma[comp]);
        div(npa+arete)-=psc;
      }
  }
#else
#endif
 
// Divide by the volume
void Op_Div_VEFP1B_Elem::volumique_P0(DoubleTab& tab_div) const
{
  const Domaine_VEF& domaine_VEF = le_dom_vef.valeur();
  int nb_elem = domaine_VEF.domaine().nb_elem_tot();
  CDoubleArrView vol = domaine_VEF.volumes().view_ro();
  DoubleArrView div = static_cast<DoubleVect&>(tab_div).view_rw();
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_elem, KOKKOS_LAMBDA(const int i) { div(i) /= vol(i); });
  end_gpu_timer(__KERNEL_NAME__);
}
 
void Op_Div_VEFP1B_Elem::volumique(DoubleTab& tab_div) const
{
  const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, le_dom_vef.valeur());
  int n = 0;
  if (domaine_VEF.get_alphaE())
    {
      volumique_P0(tab_div);
      n += domaine_VEF.nb_elem_tot();
    }
  if (domaine_VEF.get_alphaS())
    {
      int size_tot = domaine_VEF.volume_aux_sommets().size_totale();
      CDoubleArrView vol = domaine_VEF.volume_aux_sommets().view_ro();
      DoubleArrView div = static_cast<DoubleVect&>(tab_div).view_rw();
      Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), size_tot, KOKKOS_LAMBDA(const int i) { div(n + i) /= vol(i); });
      end_gpu_timer(__KERNEL_NAME__);
      n += domaine_VEF.nb_som_tot();
    }
  if (domaine_VEF.get_alphaA())
    {
      const DoubleVect& vol = domaine_VEF.get_volumes_aretes();
      int size_tot = vol.size_totale();
      ToDo_Kokkos("critical");
      for (int i = 0; i < size_tot; i++)
        tab_div(n + i) /= vol(i);
    }
}
 
void Op_Div_VEFP1B_Elem::degres_liberte() const
{
  // Set to zero the divergence at vertices with no degree of freedom
  // (vertex only shared by faces with Dirichlet boundary conditions)
  const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, le_dom_vef.valeur());
  const Domaine& domaine = domaine_VEF.domaine();
  const IntTab& elem_faces = domaine_VEF.elem_faces();
  const IntTab& face_voisins = domaine_VEF.face_voisins();
  const int nse = domaine.nb_som_elem();
  int nb_faces_tot = domaine_VEF.nb_faces_tot();
 
  Sortie& journal = Process::Journal();
  int afficher_message = 0;
  int nb_som = domaine.nb_som();
 
  int ecrire_decoupage_som = 0;
  int decoupage_som = 0;
  // Write the decoupage_som file only once
  // and only in sequential mode
  if ((Process::is_sequential()) && equation().schema_temps().nb_pas_dt() == 0)
    decoupage_som = 1;
 
  SChaine decoup_som;
  decoup_som << "1" << finl;
  decoup_som << Objet_U::dimension << " " << nb_som << finl;
  ArrOfInt somm(dimension + 2);
  for (int k = 0; k < nb_som; k++)
    {
      int sommet = domaine.get_renum_som_perio(k);
      if (nb_degres_liberte_(sommet) != 0)
        continue;
      if (!afficher_message)
        {
          afficher_message = 1;
          Cerr << finl << "Problem with the mesh used for the VEF P1Bulle discretization." << finl;
          journal << "List of nodes with no degrees of freedom :" << finl;
        }
      const double x = domaine.coord(sommet, 0);
      const double y = domaine.coord(sommet, 1);
      const double z = (Objet_U::dimension == 3) ? domaine.coord(sommet, 2) : 0.;
 
      journal << "Error node " << sommet << " ( " << x << " " << y << " " << z << " )\n";
      // Display the list of indices of real and virtual elements containing
      // this vertex. Display the letter "v" for virtual elements.
      journal << "Elements ";
      const int nb_elem_tot = domaine.nb_elem_tot();
      const int nb_elem = domaine.nb_elem();
      const IntTab& som_elem = domaine.les_elems();
      for (int elem = 0; elem < nb_elem_tot; elem++)
        for (int som = 0; som < nse; som++)
          if (som_elem(elem, som) == sommet)
            {
              journal << elem << ((elem >= nb_elem) ? "v " : " ");
 
              // Write to the decoupage_som file
              int face_opp = elem_faces(elem, som);
              int elem_opp;
              somm = -1;
              somm(0) = sommet;
 
              int elem1 = face_voisins(face_opp, 0);
              int elem2 = face_voisins(face_opp, 1);
 
              if (elem1 == elem)
                elem_opp = elem2;
              else
                elem_opp = elem1;
 
              int i = 2;
              for (int som1 = 0; som1 < nse; som1++)  // iterate over vertices of elem_opp
                {
                  int ok = 1;
                  for (int som2 = 0; som2 < nse; som2++)  // iterate over vertices of elem
                    if (som_elem(elem, som2) == som_elem(elem_opp, som1))
                      ok = 0;
                  if (ok)
                    somm(1) = som_elem(elem_opp, som1);
                  else
                    {
                      somm(i) = som_elem(elem_opp, som1);  // vertices of the shared face
                      i++;
                    }
                }
              if (decoupage_som)
                {
                  ecrire_decoupage_som = 1;
                  for (int j = 0; j < dimension + 2; j++)
                    decoup_som << somm(j) << " ";
                  decoup_som << elem << " " << elem_opp << finl;
                }
            }
      journal << "\n";
      // Display the list of faces containing this vertex.
      // For boundary faces, display the boundary condition,
      // for virtual faces, the letter "v"
      journal << "\nFaces ";
      const int nb_faces = domaine_VEF.nb_faces();
      const int nb_som_face = domaine_VEF.face_sommets().dimension(1);
      for (int face = 0; face < nb_faces_tot; face++)
        {
          for (int som = 0; som < nb_som_face; som++)
            {
              if (domaine_VEF.face_sommets(face, som) == sommet)
                {
                  journal << face;
                  if (face >= nb_faces) // Virtual face
                    journal << "v";
                  const int cl = find_cl_face(domaine, face);
                  // Real boundary face:
                  if (cl >= 0)
                    {
                      const Nom& nom_bord = domaine.frontiere(cl).le_nom();
                      journal << "(boundary=" << nom_bord << ")";
                    }
                  journal << " ";
                }
            }
        }
      journal << finl;
    }
 
  if (ecrire_decoupage_som)
    {
      decoup_som << "-1";
      SFichier fic(Objet_U::nom_du_cas() + ".decoupage_som");
      fic << decoup_som.get_str() << finl;
      fic.close();
 
    }
  if (afficher_message)
    {
      Nom nom_fichier(nom_du_cas());
      Cerr << "Look at the .log file of processor " << Process::me() << " to know which nodes" << finl;
      Cerr << "do not have enough degrees of freedom. You can also visualize" << finl;
      Cerr << "your mesh to identify which elements have a problem and" << finl;
      Cerr << "remesh your domain close to these nodes so that these have" << finl;
      Cerr << "at least a degree of freedom." << finl;
      Cerr << "The VEF P1Bulle discretization is not compatible with this type" << finl;
      Cerr << "of mesh and boundary conditions used (Dirichlet, Symmetry ,...)" << finl << finl;
      Cerr << "Other possibility, in your data file:" << finl;
      if (dimension == 2)
        Cerr << "If you have used \"Trianguler\", substituted by \"Trianguler_H\"." << finl << finl;
      if (dimension == 3)
        Cerr << "If you have used \"Tetraedriser\", substituted by \"Tetraedriser_homogene\" or \"Tetraedriser_par_prisme\"." << finl << finl;
      Cerr << "Or insert the line:" << finl;
      Cerr << "VerifierCoin " << domaine.le_nom() << " { [Read_file " << nom_fichier << ".decoupage_som] }" << finl;
      Cerr << "after the mesh is finished to be read and built." << finl;
      if (Process::is_parallel())
        Cerr << "and BEFORE the keyword \"Decouper\" during the partition of the mesh." << finl;
      else
        Cerr << "and BEFORE the keyword \"Discretiser\"." << finl;
      Cerr << "A few cells will be divided into 3 (2D) or 4 (3D)." << finl;
      Cerr << finl;
      Process::exit();
    }
}
 
int Op_Div_VEFP1B_Elem::impr(Sortie& os) const
{
  const int impr_bord = (le_dom_vef->domaine().bords_a_imprimer().est_vide() ? 0 : 1);
  const Schema_Temps_base& sch = equation().probleme().schema_temps();
  double temps = sch.temps_courant();
 
  int nb_compo = flux_bords_.dimension(1);
  // Loop over boundaries to sum the fluxes per boundary into the flux_bord array
  DoubleVect bilan(nb_compo);
  bilan = 0;
  int nb_cl = le_dom_vef->nb_front_Cl();
  // flux_bords contains the sum of flux on each boundary:
  DoubleTrav tab_flux_bords(3, nb_cl, nb_compo);
  for (int num_cl = 0; num_cl < nb_cl; num_cl++)
    {
      const Cond_lim& la_cl = la_zcl_vef->les_conditions_limites(num_cl);
      const Front_VF& frontiere_dis = ref_cast(Front_VF, la_cl->frontiere_dis());
      int ndeb = frontiere_dis.num_premiere_face();
      int nfin = ndeb + frontiere_dis.nb_faces();
      int perio = (sub_type(Periodique,la_cl.valeur()) ? 1 : 0);
      for (int face = ndeb; face < nfin; face++)
        for (int k = 0; k < nb_compo; k++)
          {
            tab_flux_bords(0, num_cl, k) += flux_bords_(face, k);
            if (perio)
              {
                if (face < (ndeb + frontiere_dis.nb_faces() / 2))
                  tab_flux_bords(1, num_cl, k) += flux_bords_(face, k);
                else
                  tab_flux_bords(2, num_cl, k) += flux_bords_(face, k);
              }
          }
    }
  // Sum on all CPUs:
  mp_sum_for_each_item(tab_flux_bords);
 
  // Print
  if (je_suis_maitre())
    {
      ouvrir_fichier(Flux_div, "", 1);
      Flux_div.add_col(temps);
      for (int num_cl = 0; num_cl < nb_cl; num_cl++)
        {
          const Cond_lim& la_cl = la_zcl_vef->les_conditions_limites(num_cl);
          int perio = (sub_type(Periodique,la_cl.valeur()) ? 1 : 0);
          for (int k = 0; k < nb_compo; k++)
            {
              if (perio)
                {
                  Flux_div.add_col(tab_flux_bords(1, num_cl, k));
                  Flux_div.add_col(tab_flux_bords(2, num_cl, k));
                }
              else
                Flux_div.add_col(tab_flux_bords(0, num_cl, k));
              bilan(k) += tab_flux_bords(0, num_cl, k);
            }
        }
 
      for (int k = 0; k < nb_compo; k++)
        Flux_div.add_col(bilan(k));
      Flux_div << finl;
    }
 
  const LIST(Nom) &Liste_bords_a_imprimer = le_dom_vef->domaine().bords_a_imprimer();
  if (!Liste_bords_a_imprimer.est_vide())
    {
      EcrFicPartage Flux_face;
      ouvrir_fichier_partage(Flux_face, "", impr_bord);
      for (int num_cl = 0; num_cl < nb_cl; num_cl++)
        {
          const Frontiere_dis_base& la_fr = la_zcl_vef->les_conditions_limites(num_cl)->frontiere_dis();
          const Cond_lim& la_cl = la_zcl_vef->les_conditions_limites(num_cl);
          const Front_VF& frontiere_dis = ref_cast(Front_VF, la_cl->frontiere_dis());
          int ndeb = frontiere_dis.num_premiere_face();
          int nfin = ndeb + frontiere_dis.nb_faces();
          if (le_dom_vef->domaine().bords_a_imprimer().contient(la_fr.le_nom()))
            {
              Flux_face << "# Flux par face sur " << la_fr.le_nom() << " au temps " << temps << " : " << finl;
              for (int face = ndeb; face < nfin; face++)
                {
                  if (dimension == 2)
                    Flux_face << "# Face a x= " << le_dom_vef->xv(face, 0) << " y= " << le_dom_vef->xv(face, 1) << " flux=";
                  else if (dimension == 3)
                    Flux_face << "# Face a x= " << le_dom_vef->xv(face, 0) << " y= " << le_dom_vef->xv(face, 1) << " z= " << le_dom_vef->xv(face, 2) << " flux=";
                  for (int k = 0; k < nb_compo; k++)
                    Flux_face << " " << flux_bords_(face, k);
                  Flux_face << finl;
                }
              Flux_face.syncfile();
            }
        }
    }
  return 1;
}