Domaine_Cl_VEF.cpp

/****************************************************************************
* Copyright (c) 2025, 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 <Dirichlet_paroi_defilante.h>
#include <Scalaire_impose_paroi.h>
#include <Check_espace_virtuel.h>
#include <Dirichlet_paroi_fixe.h>
#include <Neumann_sortie_libre.h>
#include <Domaine_VEF.h>
#include <Dirichlet_homogene.h>
#include <Navier_Stokes_std.h>
#include <MD_Vector_tools.h>
#include <Domaine_Cl_VEF.h>
#include <Champ_P0_VEF.h>
#include <Hexaedre_VEF.h>
#include <Domaine_VEF.h>
#include <Quadri_VEF.h>
#include <Champ_P1NC.h>
#include <Periodique.h>
#include <Conduction.h>
#include <Dirichlet.h>
#include <Tetra_VEF.h>
#include <Hexa_VEF.h>
#include <Symetrie.h>
#include <Domaine.h>
#include <Tri_VEF.h>
#include <Debog.h>
 
Implemente_instanciable(Domaine_Cl_VEF, "Domaine_Cl_VEF", Domaine_Cl_dis_base);
 
int Deux_Puissance(int n)
{
  switch(n)
    {
    case 0:
      return 1;
    case 1:
      return 2;
    case 2:
      return 4;
    case 3:
      return 8;
    default:
      return 8 * Deux_Puissance((int) (n - 3));
    }
}
 
Sortie& Domaine_Cl_VEF::printOn(Sortie& os) const
{
  os << volumes_entrelaces_Cl_;
  os << normales_facettes_Cl_;
  os << type_elem_Cl_;
  return os;
}
 
Entree& Domaine_Cl_VEF::readOn(Entree& is) { return Domaine_Cl_dis_base::readOn(is); }
 
void Domaine_Cl_VEF::associer(const Domaine_dis_base& dom_dis)
{
  const Domaine_VEF& le_dom_VEF = ref_cast(Domaine_VEF, dom_dis);
  int nb_faces_non_std = le_dom_VEF.nb_faces_non_std();
  const Elem_VEF_base& type_elem = le_dom_VEF.type_elem();
  int nb_fa7_elem = type_elem.nb_facette();
 
  {
    IntVect renum;
    le_dom_VEF.creer_tableau_faces(renum, RESIZE_OPTIONS::NOCOPY_NOINIT);
    renum = -1;
    // Mark real non-standard faces (first faces in the domain)
    for (int i = 0; i < nb_faces_non_std; i++)
      renum[i] = 0;
    // Exchange virtual space => mark virtual non-standard faces
    renum.echange_espace_virtuel();
    MD_Vector md_vect;
    MD_Vector_tools::creer_md_vect_renum_auto(renum, md_vect);
    MD_Vector_tools::creer_tableau_distribue(md_vect, volumes_entrelaces_Cl_, RESIZE_OPTIONS::NOCOPY_NOINIT);
  }
  {
    // Build the descriptor for the "Cl" element array:
    MD_Vector md_vect;
    MD_Vector_tools::creer_md_vect_renum(le_dom_VEF.rang_elem_non_std(), md_vect);
 
    // Create arrays for "Cl" elements:
    normales_facettes_Cl_.resize(0, nb_fa7_elem, dimension);
    MD_Vector_tools::creer_tableau_distribue(md_vect, normales_facettes_Cl_, RESIZE_OPTIONS::NOCOPY_NOINIT);
 
    vecteur_face_facette_Cl_.resize(0, nb_fa7_elem, dimension, 2);
    MD_Vector_tools::creer_tableau_distribue(md_vect, vecteur_face_facette_Cl_, RESIZE_OPTIONS::NOCOPY_NOINIT);
 
    MD_Vector_tools::creer_tableau_distribue(md_vect, type_elem_Cl_, RESIZE_OPTIONS::NOCOPY_NOINIT);
  }
}
 
/*! @brief Fill arrays.
 *
 */
void Domaine_Cl_VEF::completer(const Domaine_dis_base& un_domaine_dis)
{
  if (sub_type(Domaine_VEF, un_domaine_dis))
    {
      const Domaine_VEF& le_dom_VEF = ref_cast(Domaine_VEF, un_domaine_dis);
      remplir_type_elem_Cl(le_dom_VEF);
      remplir_volumes_entrelaces_Cl(le_dom_VEF);
      remplir_normales_facettes_Cl(le_dom_VEF);
    }
  else
    {
      Cerr << "Domaine_Cl_VEF::completer() expects a Domaine_VEF argument " << finl;
      exit();
    }
}
 
/*! @brief Called by completer(): fills type_elem_Cl_ and volumes_entrelaces_Cl_.
 *
 */
void Domaine_Cl_VEF::remplir_volumes_entrelaces_Cl(const Domaine_VEF& le_dom_VEF)
{
  Cerr << "Entering Domaine_Cl_VEF::remplir_volumes_entrelaces" << finl;
  const Domaine_VEF& z = ref_cast(Domaine_VEF, le_dom_VEF);
 
  // Extend control volumes on non-standard faces. Only in conduction for now or in P0-only mode
  if ((z.get_modif_div_face_dirichlet()) || sub_type(Conduction, equation()) || (z.get_alphaE() && !z.get_alphaS() && !z.get_alphaA()))
    {
      const DoubleVect& volumes_entrelaces = le_dom_VEF.volumes_entrelaces();
      const IntVect& rang_elem = le_dom_VEF.rang_elem_non_std();
 
      // Initialise the volumes_entrelaces_Cl array
      // non-standard faces are not a priori touched by boundary conditions
 
      for (int i = 0; i < le_dom_VEF.nb_faces_non_std(); i++)
        volumes_entrelaces_Cl_[i] = volumes_entrelaces(i);
 
      // Compute the interlaced volumes modified by boundary conditions:
 
      int nb_poly_tot = le_dom_VEF.domaine().nb_elem_tot();
      ArrOfInt poly_fait(nb_poly_tot);
      poly_fait = 0;
      for (int i = 0; i < les_conditions_limites_.size(); i++)
        {
          const Cond_lim_base& la_cl = les_conditions_limites_[i].valeur();
          if ((sub_type(Dirichlet, la_cl)) || (sub_type(Dirichlet_homogene, la_cl)))
            {
              const Front_VF& le_bord = le_dom_VEF.front_VF(i);
              const Frontiere& front = le_bord.frontiere();
              int premiere = le_bord.num_premiere_face();
              int derniere = premiere + le_bord.nb_faces();
              for (int j = premiere; j < derniere; j++)
                {
                  volumes_entrelaces_Cl_[j] = 0;
                  int elem = le_dom_VEF.face_voisins(j, 0);
                  int n_poly = rang_elem(elem);
                  if (!poly_fait[n_poly])
                    {
                      poly_fait[n_poly] = 1;
                      const Elem_VEF_base& type_elem = le_dom_VEF.type_elem();
                      type_elem.modif_volumes_entrelaces(j, elem, le_dom_VEF, volumes_entrelaces_Cl(), type_elem_Cl(n_poly));
                    }
                }
              // virtual boundary faces:
              const ArrOfInt& faces_virt = front.get_faces_virt();
              int nb_faces_virt = faces_virt.size_array();
              for (int ind_face = 0; ind_face < nb_faces_virt; ind_face++)
                {
                  int num_face = faces_virt[ind_face];
                  int elem = le_dom_VEF.face_voisins(num_face, 0);
                  int n_poly = rang_elem(elem);
                  if (n_poly != -1)
                    {
                      if (!poly_fait[n_poly])
                        {
                          poly_fait[n_poly] = 1;
                          const Elem_VEF_base& type_elem = le_dom_VEF.type_elem();
                          type_elem.modif_volumes_entrelaces_faces_joints(num_face, elem, le_dom_VEF, volumes_entrelaces_Cl(), type_elem_Cl(n_poly));
                        }
                    }
                }
            }
        }
    }
  else
    {
      const DoubleVect& volumes_entrelaces = le_dom_VEF.volumes_entrelaces();
      // Initialization of the volumes_entrelaces_Cl array: non-standard faces are a priori not affected by boundary conditions
      for (int i = 0; i < le_dom_VEF.nb_faces_non_std(); i++)
        volumes_entrelaces_Cl(i) = volumes_entrelaces(i);
      Cerr << "End of virtual space creation for volumes_entrelaces_Cl_" << finl;
    }
}
 
/*! @brief Called by completer(): fills normales_facettes_Cl_ and vecteur_face_facette_Cl_.
 *
 *  RENAME THIS METHOD
 *
 */
void Domaine_Cl_VEF::remplir_normales_facettes_Cl(const Domaine_VEF& le_dom_VEF)
{
  const Domaine& z = le_dom_VEF.domaine();
  const Domaine& dom = z;
  const DoubleTab& les_coords = dom.coord_sommets();
  const IntTab& les_Polys = z.les_elems();
  const Elem_VEF_base& elemvef = le_dom_VEF.type_elem();
  const IntVect& rang_elem = le_dom_VEF.rang_elem_non_std();
  const IntTab& elem_faces = le_dom_VEF.elem_faces();
  const DoubleTab& xv = le_dom_VEF.xv();
  const IntTab& KEL = elemvef.KEL();
 
  int num_elem, fa7;
  int isom, ncomp;
  int idirichlet = -1, n1 = -1, n2 = -1, n3 = -1;
  int nb_elem = z.nb_elem();
  int nsom = z.nb_som_elem();
  int nfa7 = elemvef.nb_facette();
  //int nb_som_facette = domaine().type_elem().nb_som_face();
  //if ( sub_type(Hexaedre_VEF,domaine().type_elem().valeur())) { nb_som_facette--; }
  int nb_som_facette = dimension;
 
  IntVect num_som(nsom);
  DoubleTab x(nsom, dimension);
  DoubleVect xg(dimension);
  //DoubleVect xj0(dimension);
  //DoubleVect u(dimension);
  //DoubleVect v(dimension);
 
  // Compute the facet normals modified by boundary conditions:
 
  for (int elem = 0; elem < nb_elem; elem++)
    {
 
      if (rang_elem(elem) != -1)
        {
          num_elem = rang_elem(elem);
 
          // index of the vertices of the polyhedron:
          for (isom = 0; isom < nsom; isom++)
            num_som[isom] = les_Polys(elem, isom);
 
          // Compute the coordinates of the polyhedron vertices:
          for (isom = 0; isom < nsom; isom++)
            for (ncomp = 0; ncomp < dimension; ncomp++)
              x(isom, ncomp) = les_coords(num_som[isom], ncomp);
 
          // Compute the coordinates of point G based on the polyhedron type
          // idirichlet = number of Dirichlet faces of the element
          // for triangle if idirichlet=2, n1=vertex coinciding with G
          // for quadrilateral refer to Quadri_VEF.cpp
 
          elemvef.calcul_xg(xg, x, type_elem_Cl_[num_elem], idirichlet, n1, n2, n3);
          // Compute values of the normales_facettes_Cl array:
 
          for (fa7 = 0; fa7 < nfa7; fa7++)
            {
              elemvef.creer_normales_facettes_Cl(normales_facettes_Cl(), fa7, num_elem, x, xg, z);
              // Correction of the normales_facettes_Cl array values:
              elemvef.modif_normales_facettes_Cl(normales_facettes_Cl(), fa7, num_elem, idirichlet, n1, n2, n3);
              int num1 = elem_faces(elem, KEL(0, fa7));
              int num2 = elem_faces(elem, KEL(1, fa7));
              for (int dir = 0; dir < dimension; dir++)
                {
                  double coord_centre_fa7 = xg(dir);
 
                  for (int num_som_fa7 = 0; num_som_fa7 < nb_som_facette - 1; num_som_fa7++)
                    {
                      int isom_loc = KEL(num_som_fa7 + 2, fa7);
                      int isom_glob = les_Polys(elem, isom_loc);
                      coord_centre_fa7 += les_coords(isom_glob, dir);
                    }
                  coord_centre_fa7 /= nb_som_facette;
 
                  vecteur_face_facette_Cl_(num_elem, fa7, dir, 0) = coord_centre_fa7 - xv(num1, dir);
                  vecteur_face_facette_Cl_(num_elem, fa7, dir, 1) = coord_centre_fa7 - xv(num2, dir);
                }
 
            }
        }
    }
  normales_facettes_Cl_.echange_espace_virtuel();
  vecteur_face_facette_Cl_.echange_espace_virtuel();
}
 
int trois_puissance(int n)
{
  switch(n)
    {
    case 0:
      return 1;
    case 1:
      return 3;
    case 2:
      return 9;
    case 3:
      return 27;
    default:
      return 27 * trois_puissance(n - 3);
    }
}
 
/*! @brief Called by remplir_volumes_entrelaces_Cl(): fills type_elem_Cl_.
 *
 */
void Domaine_Cl_VEF::remplir_type_elem_Cl(const Domaine_VEF& le_dom_VEF)
{
  const Domaine& z = le_dom_VEF.domaine();
  int nfac = z.type_elem()->nb_faces();   // in Elem_geom
  const IntTab& elem_faces = le_dom_VEF.elem_faces();
  const IntTab& face_voisins = le_dom_VEF.face_voisins();
  const IntVect& rang_elem = le_dom_VEF.rang_elem_non_std();
  type_elem_Cl_ = 0;
  int elem, num_elem;
 
  for (int i = 0; i < les_conditions_limites_.size(); i++)
    {
      Cond_lim_base& la_cl = les_conditions_limites_[i].valeur();
      if ((sub_type(Dirichlet, la_cl)) || (sub_type(Dirichlet_homogene, la_cl)))
        {
          const Front_VF& le_bord = le_dom_VEF.front_VF(i);
          int premiere = le_bord.num_premiere_face();
          int derniere = premiere + le_bord.nb_faces();
          for (int j = premiere; j < derniere; j++)
            {
              elem = face_voisins(j, 0);
              int numero1, numero2 = -1;
              for (numero1 = 0; numero1 < nfac; numero1++)
                if (elem_faces(elem, numero1) == j)
                  numero2 = numero1;
              if (numero2 == -1)
                Cerr << "Domaine_Cl_VEF::creer_type_elem_Cl() NOT POSSIBLE!! " << finl;
              num_elem = rang_elem(elem);
              assert(num_elem != -1);
              if (sub_type(Tri_VEF,le_dom_VEF.type_elem()) || sub_type(Tetra_VEF, le_dom_VEF.type_elem()))
                type_elem_Cl_[num_elem] += Deux_Puissance((int) (nfac - 1 - numero2));
              else if (sub_type(Quadri_VEF,le_dom_VEF.type_elem()) || sub_type(Hexa_VEF, le_dom_VEF.type_elem()))
                type_elem_Cl_[num_elem] += trois_puissance((int) (nfac - 1 - numero2));
              else
                {
                  Cerr << "error in Domaine_Cl_VEF::remplir_type_elem_Cl wrong element type " << finl;
                  exit();
                }
 
            }
        }
    }
 
  type_elem_Cl_.echange_espace_virtuel();
}
 
/*! @brief Imposes boundary conditions at the time "temps" of the Champ_Inc.
 *
 */
void Domaine_Cl_VEF::imposer_cond_lim(Champ_Inc_base& ch, double temps)
{
  DoubleTab& ch_tab = ch.valeurs(temps);
  int nb_comp = ch.nb_comp();
  if (sub_type(Champ_P0_VEF, ch)) { /* Do nothing */ }
  else if (sub_type(Champ_P1NC,ch) || sub_type(Champ_Q1NC, ch))
    {
      for (int i = 0; i < nb_cond_lim(); i++)
        {
          const Cond_lim_base& la_cl = les_conditions_limites(i).valeur();
          const Front_VF& le_bord = ref_cast(Front_VF, la_cl.frontiere_dis());
          int ndeb = le_bord.num_premiere_face();
          int nfin = ndeb + le_bord.nb_faces();
          if (sub_type(Periodique, la_cl))
            {
              // Ensure the field has the same value
              // on two periodic faces that are opposite each other
              const Periodique& la_cl_perio = ref_cast(Periodique, la_cl);
              CIntArrView face_associee = la_cl_perio.face_associee().view_ro();
              if (nb_comp==1)
                {
                  DoubleArrView tab = static_cast<ArrOfDouble&>(ch_tab).view_wo();
                  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin), KOKKOS_LAMBDA(const int num_face)
                  {
                    int voisine = face_associee(num_face - ndeb) + ndeb;
                    if (tab(num_face) != tab(voisine))
                      {
                        double moy = 0.5 * (tab(num_face) + tab(voisine));
                        tab(num_face) = moy;
                        tab(voisine) = moy;
                      }
                  });
                }
              else
                {
                  DoubleTabView tab = ch_tab.view_wo();
                  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin), KOKKOS_LAMBDA(const int num_face)
                  {
                    int voisine = face_associee(num_face - ndeb) + ndeb;
                    for (int ncomp = 0; ncomp < nb_comp; ncomp++)
                      if (tab(num_face, ncomp) != tab(voisine, ncomp))
                        {
                          double moy = 0.5 * (tab(num_face, ncomp) + tab(voisine, ncomp));
                          tab(num_face, ncomp) = moy;
                          tab(voisine, ncomp) = moy;
                        }
                  });
                }
              end_gpu_timer(__KERNEL_NAME__);
            }
          else if ((sub_type(Symetrie, la_cl)) && (ch.nature_du_champ() == vectoriel))
            {
              if (nb_comp == dimension)
                {
                  CDoubleTabView face_normales = domaine_vef().face_normales().view_ro();
                  DoubleTabView tab = ch_tab.view_rw();
                  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin), KOKKOS_LAMBDA(const int num_face)
                  {
                    double surf = 0, flux = 0;
                    for (int ncomp = 0; ncomp < nb_comp; ncomp++)
                      {
                        double face_n = face_normales(num_face, ncomp);
                        flux += tab(num_face, ncomp) * face_n;
                        surf += face_n * face_n;
                      }
                    // flux /= surf; // Fixed bug: Arithmetic exception
                    if (std::fabs(surf) >= DMINFLOAT)
                      flux /= surf;
                    for (int ncomp = 0; ncomp < nb_comp; ncomp++)
                      tab(num_face, ncomp) -= flux * face_normales(num_face, ncomp);
                  });
                  end_gpu_timer(__KERNEL_NAME__);
                }
            }
          else if (sub_type(Dirichlet_entree_fluide, la_cl))
            {
              const Dirichlet_entree_fluide& la_cl_diri = ref_cast(Dirichlet_entree_fluide, la_cl);
              CDoubleTabView val_imp = la_cl_diri.tab_val_imp(temps).view_ro();
              if (nb_comp == 1)
                {
                  DoubleArrView tab = static_cast<ArrOfDouble&>(ch_tab).view_wo();
                  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin),
                                       KOKKOS_LAMBDA(const int num_face)
                  {
                    tab(num_face) = val_imp(num_face - ndeb, 0);
                  });
                }
              else
                {
                  DoubleTabView tab = ch_tab.view_wo();
                  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin), KOKKOS_LAMBDA(const int num_face)
                  {
                    for (int ncomp = 0; ncomp < nb_comp; ncomp++)
                      tab(num_face, ncomp) = val_imp(num_face - ndeb, ncomp);
                  });
                }
              end_gpu_timer(__KERNEL_NAME__);
            }
          else if (sub_type(Scalaire_impose_paroi, la_cl))
            {
              const Scalaire_impose_paroi& la_cl_diri = ref_cast(Scalaire_impose_paroi, la_cl);
              CDoubleTabView val_imp = la_cl_diri.tab_val_imp(temps).view_ro();
              if (nb_comp == 1)
                {
                  DoubleArrView tab = static_cast<ArrOfDouble&>(ch_tab).view_wo();
                  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin),
                                       KOKKOS_LAMBDA(const int num_face)
                  {
                    tab(num_face) = val_imp(num_face - ndeb, 0);
                  });
                }
              else
                {
                  DoubleTabView tab = ch_tab.view_wo();
                  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin), KOKKOS_LAMBDA(const int num_face)
                  {
                    for (int ncomp = 0; ncomp < nb_comp; ncomp++)
                      tab(num_face, ncomp) = val_imp(num_face - ndeb, ncomp);
                  });
                }
              end_gpu_timer(__KERNEL_NAME__);
            }
          else if ((sub_type(Dirichlet_paroi_fixe, la_cl)) && (ch.nature_du_champ() == multi_scalaire))
            {
              if (nb_comp == 1)
                {
                  DoubleArrView tab = static_cast<ArrOfDouble&>(ch_tab).view_wo();
                  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin),
                                       KOKKOS_LAMBDA(const int num_face)
                  {
                    tab(num_face) = 0;
                  });
                }
              else
                {
                  DoubleTabView tab = ch_tab.view_wo();
                  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin), KOKKOS_LAMBDA(const int num_face)
                  {
                    for (int ncomp = 0; ncomp < nb_comp; ncomp++)
                      {
                        tab(num_face, 0) = 0;
                        tab(num_face, 1) = 0;
                      }
                  });
                }
              end_gpu_timer(__KERNEL_NAME__);
            }
          else if (sub_type(Dirichlet_paroi_fixe, la_cl))
            {
              if (nb_comp == 1)
                {
                  DoubleArrView tab = static_cast<ArrOfDouble&>(ch_tab).view_wo();
                  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin),
                                       KOKKOS_LAMBDA(const int num_face)
                  {
                    tab(num_face) = 0;
                  });
                }
              else
                {
                  DoubleTabView tab = ch_tab.view_wo();
                  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin), KOKKOS_LAMBDA(const int num_face)
                  {
                    for (int ncomp = 0; ncomp < nb_comp; ncomp++)
                      tab(num_face, ncomp) = 0;
                  });
                }
              end_gpu_timer(__KERNEL_NAME__);
            }
          else if (sub_type(Dirichlet_paroi_defilante, la_cl))
            {
              const Dirichlet_paroi_defilante& la_cl_diri = ref_cast(Dirichlet_paroi_defilante, la_cl);
              CDoubleTabView face_normales = domaine_vef().face_normales().view_ro();
              if (nb_comp != dimension)
                {
                  Cerr << "Case not handled in Domaine_Cl_VEF::imposer_cond_lim." << finl;
                  exit();
                }
              else
                {
                  CDoubleTabView val_imp = la_cl_diri.tab_val_imp(temps).view_ro();
                  DoubleTabView tab = ch_tab.view_wo();
                  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin), KOKKOS_LAMBDA(const int num_face)
                  {
                    double surf = 0, flux = 0;
                    for (int ncomp = 0; ncomp < nb_comp; ncomp++)
                      {
                        double face_n = face_normales(num_face, ncomp);
                        flux += val_imp(num_face - ndeb, ncomp) * face_n;
                        surf += face_n * face_n;
                      }
                    // flux /= surf; // Fixed bug: Arithmetic exception
                    if (std::fabs(surf) >= DMINFLOAT)
                      flux /= surf;
                    for (int ncomp = 0; ncomp < nb_comp; ncomp++)
                      tab(num_face, ncomp) = val_imp(num_face - ndeb, ncomp) -
                                             flux * face_normales(num_face, ncomp);
                  });
                  end_gpu_timer(__KERNEL_NAME__);
                }
            }
        }
    }
  else
    {
      Cerr << "The type of OWN_PTR(Champ_Inc_base) " << ch.que_suis_je() << " is not handled in VEF\n" << finl;
      exit();
    }
  ch_tab.echange_espace_virtuel();
 
  // PRESSURE PART
  // In the Navier-Stokes case and if the pressure BC is imposed strongly at vertices, impose the pressure values at vertices
  if (sub_type(Navier_Stokes_std, ch.equation()))
    {
      const Domaine_VEF& domaine_vef = ref_cast(Domaine_VEF, equation().domaine_dis());
      if ((domaine_vef.get_cl_pression_sommet_faible() == 0) && (domaine_vef.get_alphaS()))
        {
          int nps = domaine_vef.numero_premier_sommet();
          int nbsom_tot = domaine_vef.nb_som_tot();
          DoubleTrav tab_surf_loc(nbsom_tot);
          tab_surf_loc = 0;
          DoubleTab& tab_pression = ref_cast(Navier_Stokes_std,ch.equation()).pression().valeurs();
          int nb_cond_lims = nb_cond_lim();
          // First loop: zero out the pressure at vertices
          for (int i = 0; i < nb_cond_lims; i++)
            {
              const Cond_lim_base& la_cl = les_conditions_limites(i).valeur();
              if (sub_type(Neumann_sortie_libre, la_cl))
                {
                  const Front_VF& le_bord = ref_cast(Front_VF, la_cl.frontiere_dis());
                  int nbsf = domaine_vef.face_sommets().dimension(1);
                  int num2 = le_bord.nb_faces_tot();
                  CIntTabView faces = domaine_vef.face_sommets().view_ro();
                  CIntArrView num_face = le_bord.num_face().view_ro();
                  DoubleArrView pression = static_cast<DoubleVect&>(tab_pression).view_rw();
                  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                                       Kokkos::RangePolicy<>(0, num2), KOKKOS_LAMBDA(
                                         const int ind_face)
                  {
                    int face = num_face(ind_face);
                    for (int som = 0; som < nbsf; som++)
                      {
                        int som_glob = faces(face, som);
                        Kokkos::atomic_store(&pression(nps + som_glob), 0);
                      }
                  });
                  end_gpu_timer(__KERNEL_NAME__);
                }
            }
 
          // Second loop: add the contribution of each face
          for (int i = 0; i < nb_cond_lims; i++)
            {
              const Cond_lim_base& la_cl = les_conditions_limites(i).valeur();
              if (sub_type(Neumann_sortie_libre, la_cl))
                {
                  const Neumann_sortie_libre& la_sortie_libre = ref_cast(Neumann_sortie_libre, la_cl);
                  const Front_VF& le_bord = ref_cast(Front_VF, la_cl.frontiere_dis());
                  int nbsf = domaine_vef.face_sommets().dimension(1);
                  int num2 = le_bord.nb_faces_tot();
                  CIntArrView num_face = le_bord.num_face().view_ro();
                  CIntTabView faces = domaine_vef.face_sommets().view_ro();
                  CDoubleArrView flux_impose = static_cast<const DoubleVect&>(la_sortie_libre.flux_impose(true)).view_ro();
                  CDoubleArrView face_surfaces = domaine_vef.face_surfaces().view_ro();
                  DoubleArrView surf_loc = static_cast<DoubleVect&>(tab_surf_loc).view_rw();
                  DoubleArrView pression = static_cast<DoubleVect&>(tab_pression).view_rw();
                  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                                       Kokkos::RangePolicy<>(0, num2), KOKKOS_LAMBDA(
                                         const int ind_face)
                  {
                    int face = num_face(ind_face);
                    double P_imp = flux_impose(ind_face);
                    double face_surf = face_surfaces(face);
                    for (int som = 0; som < nbsf; som++)
                      {
                        int som_glob = faces(face, som);
                        Kokkos::atomic_add(&pression(nps + som_glob), P_imp * face_surf);
                        Kokkos::atomic_add(&surf_loc[som_glob], face_surf);
                      }
                  });
                  end_gpu_timer(__KERNEL_NAME__);
                }
            }
          // Third loop: divide by the surface area
          for (int i = 0; i < nb_cond_lims; i++)
            {
              const Cond_lim_base& la_cl = les_conditions_limites(i).valeur();
              if (sub_type(Neumann_sortie_libre, la_cl))
                {
                  DoubleArrView surf_loc = static_cast<DoubleVect&>(tab_surf_loc).view_rw();
                  DoubleArrView pression = static_cast<DoubleVect&>(tab_pression).view_rw();
                  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                                       Kokkos::RangePolicy<>(0, nbsom_tot), KOKKOS_LAMBDA(
                                         const int som_glob)
                  {
                    if (surf_loc[som_glob] > 0)
                      {
                        pression(nps + som_glob) /= surf_loc[som_glob];
                        surf_loc[som_glob] = 0;
                      }
                  });
                  end_gpu_timer(__KERNEL_NAME__);
                }
            }
          tab_pression.echange_espace_virtuel();
          Debog::verifier("pression dans Domaine_Cl_VEF::imposer_cond_lim", tab_pression);
        }
    }
}
 
int Domaine_Cl_VEF::nb_faces_sortie_libre() const
{
  int compteur = 0;
  for (const auto &itr : les_conditions_limites_)
    {
      if (sub_type(Neumann_sortie_libre, itr.valeur()))
        {
          const Front_VF& le_bord = ref_cast(Front_VF, itr->frontiere_dis());
          compteur += le_bord.nb_faces();
        }
    }
  return compteur;
}
 
int Domaine_Cl_VEF::nb_bord_periodicite() const
{
  int compteur = 0;
  for (const auto &itr : les_conditions_limites_)
    {
      if (sub_type(Periodique, itr.valeur()))
        compteur++;
    }
  return compteur;
}
 
int Domaine_Cl_VEF::initialiser(double temps)
{
  Domaine_Cl_dis_base::initialiser(temps);
 
  if (nb_bord_periodicite() > 0)
    {
      if (modif_perio_fait_ == 0)
        {
          Cerr << "Updating Domaine_Cl_VEF for periodicity" << finl;
          const Domaine_VEF& le_dom_VEF = ref_cast(Domaine_VEF, domaine_dis());
          const DoubleVect& volumes_entrelaces = le_dom_VEF.volumes_entrelaces();
          remplir_volumes_entrelaces_Cl(le_dom_VEF);
          for (int i = 0; i < les_conditions_limites_.size(); i++)
            {
              const Cond_lim_base& la_cl = les_conditions_limites_[i].valeur();
 
              if (sub_type(Periodique, la_cl))
                {
 
                  const Periodique& la_cl_perio = ref_cast(Periodique, la_cl);
                  const Front_VF& le_bord = le_dom_VEF.front_VF(i);
                  int premiere = le_bord.num_premiere_face();
                  int derniere = premiere + le_bord.nb_faces();
                  IntVect fait(le_bord.nb_faces());
                  fait = 0;
                  int ind_fac_loc;
 
                  for (int j = premiere; j < derniere; j++)
                    {
                      ind_fac_loc = j - premiere;
 
                      if (fait[ind_fac_loc] == 0)
                        {
                          int fac_associee = la_cl_perio.face_associee(j - premiere) + premiere;
                          if (volumes_entrelaces_Cl_[j] != volumes_entrelaces[j])
                            volumes_entrelaces_Cl_[fac_associee] = volumes_entrelaces_Cl_[j];
                          else
                            volumes_entrelaces_Cl_[j] = volumes_entrelaces_Cl_[fac_associee];
 
                          fait[ind_fac_loc] = 1;
                          fait[fac_associee - premiere] = 1;
                        }
                    }
                }
            }
 
          modif_perio_fait_ = 1;
        }
    }
  return 1;
}
 
Domaine_VEF& Domaine_Cl_VEF::domaine_vef()
{
  return ref_cast(Domaine_VEF, domaine_dis());
}
 
const Domaine_VEF& Domaine_Cl_VEF::domaine_vef() const
{
  return ref_cast(Domaine_VEF, domaine_dis());
}