Assembleur_P_VEF.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 <Assembleur_P_VEF.h>
#include <Domaine_Cl_VEF.h>
#include <Domaine_VEF.h>
#include <Periodique.h>
#include <Neumann_sortie_libre.h>
#include <Matrice_Bloc.h>
#include <Milieu_base.h>
#include <Robin_VEF.h>
 
 
Implemente_instanciable(Assembleur_P_VEF,"Assembleur_P_VEF",Assembleur_base);
 
Sortie& Assembleur_P_VEF::printOn(Sortie& s ) const
{
  return s << que_suis_je() << " " << le_nom() ;
}
 
Entree& Assembleur_P_VEF::readOn(Entree& s )
{
  return Assembleur_base::readOn(s);
}
 
int Assembleur_P_VEF::assembler(Matrice& la_matrice)
{
  // Si rho est constant, on resout avec la pression P*=P/rho
  const DoubleVect& volumes_entrelaces_ref=le_dom_VEF->volumes_entrelaces();
  DoubleVect tab_volumes_entrelaces(volumes_entrelaces_ref);
  const DoubleVect& tab_volumes_entrelaces_cl=le_dom_Cl_VEF->volumes_entrelaces_Cl();
  int size=tab_volumes_entrelaces_cl.size();
  {
    CDoubleArrView volumes_entrelaces_cl = static_cast<const ArrOfDouble&>(tab_volumes_entrelaces_cl).view_ro();
    DoubleArrView volumes_entrelaces = static_cast<ArrOfDouble&>(tab_volumes_entrelaces).view_rw();
    Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(0, size), KOKKOS_LAMBDA(const int f)
    {
      if (volumes_entrelaces_cl(f)!=0)
        volumes_entrelaces(f)=volumes_entrelaces_cl(f);
    });
    end_gpu_timer(__KERNEL_NAME__);
  }
 
  tab_volumes_entrelaces.echange_espace_virtuel();
  // On assemble la matrice
  return assembler_mat(la_matrice,tab_volumes_entrelaces,1,1);
}
 
 
 
void calculer_inv_volume_special(DoubleTab& tab_inv_volumes_entrelaces, const Domaine_Cl_VEF& domaine_Cl_VEF,const DoubleTab& tab_volumes_entrelaces)
{
  tab_inv_volumes_entrelaces=tab_volumes_entrelaces;
  int taille=tab_volumes_entrelaces.dimension_tot(0);
  {
    CDoubleTabView volumes_entrelaces = tab_volumes_entrelaces.view_ro();
    DoubleTabView inv_volumes_entrelaces = tab_inv_volumes_entrelaces.view_rw();
    const int dimension = Objet_U::dimension;
    Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(0, taille), KOKKOS_LAMBDA(const int i)
    {
      for (int comp=0; comp<dimension; comp++)
        inv_volumes_entrelaces(i,comp)=1./volumes_entrelaces(i,comp);
    });
    end_gpu_timer(__KERNEL_NAME__);
  }
}
void Assembleur_P_VEF::calculer_inv_volume(DoubleTab& inv_volumes_entrelaces, const Domaine_Cl_VEF& domaine_Cl_VEF,const DoubleVect& volumes_entrelaces)
{
  // maintenant l 'inverse du volume est un DoubleTab
  // c'est pour faire fonctionner le Piso
  const DoubleTab* doubleT = dynamic_cast<const DoubleTab*>(&volumes_entrelaces);
  if (doubleT)
    {
      calculer_inv_volume_special(inv_volumes_entrelaces, domaine_Cl_VEF,*doubleT);
      return;
    }
  int taille=volumes_entrelaces.size_totale();
  inv_volumes_entrelaces.resize(taille,Objet_U::dimension);
  if (0)
    {
      // cette facon de calculer le volume entrelace pourra etre bonne dans
      // l'avenir a condition de reflechir au pb des porosites
      // et SURTOUT au pb du simpler ou l'on ne veut pas passer par la
      DoubleTab tmp;
      tmp=(inv_volumes_entrelaces);
      tmp=1;
      domaine_Cl_VEF.equation().solv_masse().appliquer(tmp);
      int sz=inv_volumes_entrelaces.size_totale();
      for (int i=0; i<sz; i++)
        {
          inv_volumes_entrelaces(i)=tmp(i);
          if (!est_egal(inv_volumes_entrelaces(i),1./volumes_entrelaces(i)))
            Cerr<<i<<" "<<inv_volumes_entrelaces(i)-1./volumes_entrelaces(i)<<" "<<inv_volumes_entrelaces(i)<<finl;;
        }
      Process::exit();
    }
  else
    {
      CDoubleArrView porosite_face = static_cast<const ArrOfDouble&>(equation().milieu().porosite_face()).view_ro();
      CDoubleArrView volumes_entrelaces_v = static_cast<const ArrOfDouble&>(volumes_entrelaces).view_ro();
      DoubleTabView inv_volumes_entrelaces_v = inv_volumes_entrelaces.view_rw();
      const int dim = Objet_U::dimension;
      Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(0, taille), KOKKOS_LAMBDA(const int i)
      {
        for (int comp=0; comp<dim; comp++)
          inv_volumes_entrelaces_v(i,comp)=1./volumes_entrelaces_v(i)*porosite_face(i);
      });
      end_gpu_timer(__KERNEL_NAME__);
    }
}
 
 
int Assembleur_P_VEF::assembler_mat(Matrice& la_matrice, const DoubleVect& volumes_entrelaces, int incr_pression, int resoudre_en_u)
{
  // On fixe les drapeaux de Assembleur_base
  set_resoudre_increment_pression(incr_pression);
  set_resoudre_en_u(resoudre_en_u);
  const Domaine_Cl_VEF& domaine_Cl_VEF = le_dom_Cl_VEF.valeur();
  DoubleTab inverse_quantitee_entrelacee;
  calculer_inv_volume(inverse_quantitee_entrelacee, domaine_Cl_VEF, volumes_entrelaces);
  remplir(la_matrice, inverse_quantitee_entrelacee);
  modifier_matrice(la_matrice);
  return 1;
}
 
int Assembleur_P_VEF::remplir(Matrice& la_matrice, const DoubleTab& inverse_quantitee_entrelacee)
{
  has_P_ref=0;
  // Matrice de pression :matrice creuse de taille nb_poly x nb_poly
  // Cette fonction range la matrice dans une structure de matrice morse
  // bien adaptee aux matrices creuses.
  // On commence par calculer les tailles des tableaux tab1 et tab2
  // (coeff_ a la meme taille que tab2)
  //   A chaque polyedre on associe :
  //   - une liste d'entiers voisins[i] = {j>i t.q Mij est non nul }
  //   - une liste de reels  valeurs[i] = {Mij pour j dans Voisins[i]}
  //   - un reel terme_diag
  // Implementation temporaire:
  // On assemble une matrice de pression pour une equation d'hydraulique
  // On injecte dans cette matrice les conditions aux limites
  // On peut faire cela car a priori la matrice de pression n'est pas
  // partagee par plusieurs equations sur une meme domaine.
 
  const Domaine_VEF& le_dom = le_dom_VEF.valeur();
  const Domaine_Cl_VEF& le_dom_cl = le_dom_Cl_VEF.valeur();
  les_coeff_pression.resize(le_dom_cl.nb_faces_Cl());
  int n1 = le_dom.domaine().nb_elem_tot();
  int n2 = le_dom.domaine().nb_elem();
 
 
  // Rajout des porosites.
 
  la_matrice.typer("Matrice_Bloc"); // En fait Matrice_Bloc_Sym ?
  Matrice_Bloc& matrice=ref_cast(Matrice_Bloc, la_matrice.valeur());
  matrice.dimensionner(2,2);
  matrice.get_bloc(0,0).typer("Matrice_Morse_Sym");
  matrice.get_bloc(0,1).typer("Matrice_Morse");
  matrice.get_bloc(1,0).typer("Matrice_Morse");
  matrice.get_bloc(1,1).typer("Matrice_Morse"); // En fait Matrice_Morse_Sym ?
 
  Matrice_Morse_Sym& MBrr = ref_cast(Matrice_Morse_Sym,matrice.get_bloc(0,0).valeur());
  Matrice_Morse& MBrv = ref_cast (Matrice_Morse,matrice.get_bloc(0,1).valeur());
  Matrice_Morse& MBvr = ref_cast (Matrice_Morse, matrice.get_bloc(1,0).valeur());
  Matrice_Morse& MBvv = ref_cast (Matrice_Morse, matrice.get_bloc(1,1).valeur());
 
  MBrr.dimensionner(n2,0);
  MBrv.dimensionner(n2,0);
  MBvv.dimensionner(n1-n2,0);
  // Le sous blocs vr est dimensionne et nul
  MBvr.dimensionner(n1-n2,n2,0);
  MBvr.get_set_tab1() = 1;
 
  // On traite les faces internes:
  int ndeb = le_dom_VEF->premiere_face_int();
  int nfin = le_dom_VEF->nb_faces_tot();
  int nb_faces = le_dom_VEF->nb_faces();
 
#ifdef TRUST_USE_GPU
  ArrOfTID rang_voisinRR(n2);
  ArrOfTID rang_voisinRV(n2);
  ArrOfTID rang_voisinVV(n1-n2);
#else
  ArrOfInt rang_voisinRR(n2);
  ArrOfInt rang_voisinRV(n2);
  ArrOfInt rang_voisinVV(n1-n2);
#endif
  rang_voisinRR=1; // Diagonale
  rang_voisinRV=0; // Pas de diagonale
  rang_voisinVV=1; // Diagonale
 
  CIntTabView face_voisins = le_dom.face_voisins().view_ro();
  CIntArrView ind_faces_virt_bord = le_dom_VEF->ind_faces_virt_bord().view_ro();
  auto rang_voisinRR_v = rang_voisinRR.view_rw();
  auto rang_voisinRV_v = rang_voisinRV.view_rw();
  auto rang_voisinVV_v = rang_voisinVV.view_rw();
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(ndeb, nfin), KOKKOS_LAMBDA(const int num_face)
  {
    int elem1 = face_voisins(num_face, 0);
    int elem2 = face_voisins(num_face, 1);
    const bool is_face_virt_bord = (num_face >= nb_faces) && (ind_faces_virt_bord(num_face - nb_faces) != -1);
    if (!is_face_virt_bord && elem1 != -1 && elem2 != -1)
      {
        if (elem1 > elem2)
          {
            if(elem1 < n2)
              Kokkos::atomic_add(&rang_voisinRR_v(elem2), 1);
            else
              {
                if(elem2 < n2)
                  Kokkos::atomic_add(&rang_voisinRV_v(elem2), 1);
                else
                  Kokkos::atomic_add(&rang_voisinVV_v(elem2 - n2), 1);
              }
          }
        else // elem2 >= elem1
          {
            if(elem2 < n2)
              Kokkos::atomic_add(&rang_voisinRR_v(elem1), 1);
            else
              {
                if(elem1 < n2)
                  Kokkos::atomic_add(&rang_voisinRV_v(elem1), 1);
                else
                  Kokkos::atomic_add(&rang_voisinVV_v(elem1 - n2), 1);
              }
          }
      }
  });
  end_gpu_timer(__KERNEL_NAME__);
 
  // Prise en compte des conditions de type periodicite
  const Conds_lim& les_cl = le_dom_cl.les_conditions_limites();
  for (int i=0; i<les_cl.size(); i++)
    {
      const Cond_lim& la_cl = les_cl[i];
 
      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_tot = le_bord.nb_faces_tot();
          CIntArrView front_num_face = le_bord.num_face().view_ro();
          CIntArrView face_associee = la_cl_perio.face_associee().view_ro();
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(0, nb_faces_bord_tot), KOKKOS_LAMBDA(const int ind_face)
          {
            if (ind_face < face_associee(ind_face)) // Process each periodic pair once
              {
                int num_face = front_num_face(ind_face);
                int elem1 = face_voisins(num_face, 0);
                int elem2 = face_voisins(num_face, 1);
                if (elem1 != -1 && elem2 != -1)
                  {
                    if (elem1 > elem2)
                      {
                        if(elem1 < n2)
                          Kokkos::atomic_add(&rang_voisinRR_v(elem2), 1);
                        else
                          {
                            if(elem2 < n2)
                              Kokkos::atomic_add(&rang_voisinRV_v(elem2), 1);
                            else
                              Kokkos::atomic_add(&rang_voisinVV_v(elem2 - n2), 1);
                          }
                      }
                    else // elem2 >= elem1
                      {
                        if(elem2 < n2)
                          Kokkos::atomic_add(&rang_voisinRR_v(elem1), 1);
                        else
                          {
                            if(elem1 < n2)
                              Kokkos::atomic_add(&rang_voisinRV_v(elem1), 1);
                            else
                              Kokkos::atomic_add(&rang_voisinVV_v(elem1 - n2), 1);
                          }
                      }
                  }
              }
          });
          end_gpu_timer(__KERNEL_NAME__);
        }
    }
 
  auto& tab1RR = MBrr.get_set_tab1();
  auto& tab2RR = MBrr.get_set_tab2();
  auto& tab1RV = MBrv.get_set_tab1();
  auto& tab2RV = MBrv.get_set_tab2();
  auto& tab1VV = MBvv.get_set_tab1();
  auto& tab2VV = MBvv.get_set_tab2();
 
  tab1RR(0)=1;
  tab1RV(0)=1;
  tab1VV(0)=1;
  auto tab1RR_v = tab1RR.view_rw();
  auto tab1RV_v = tab1RV.view_rw();
  auto tab1VV_v = tab1VV.view_rw();
  using tab1_value_t = typename decltype(tab1RR_v)::value_type;
  Kokkos::parallel_scan(start_gpu_timer(__KERNEL_NAME__), range_1D(0, n2), KOKKOS_LAMBDA(const int i, tab1_value_t& update, const bool final)
  {
    update += rang_voisinRR_v(i);
    if (final) tab1RR_v(i+1) = update + 1;
  });
  end_gpu_timer(__KERNEL_NAME__);
  Kokkos::parallel_scan(start_gpu_timer(__KERNEL_NAME__), range_1D(0, n2), KOKKOS_LAMBDA(const int i, tab1_value_t& update, const bool final)
  {
    update += rang_voisinRV_v(i);
    if (final) tab1RV_v(i+1) = update + 1;
  });
  end_gpu_timer(__KERNEL_NAME__);
  Kokkos::parallel_scan(start_gpu_timer(__KERNEL_NAME__), range_1D(0, n1-n2), KOKKOS_LAMBDA(const int i, tab1_value_t& update, const bool final)
  {
    update += rang_voisinVV_v(i);
    if (final) tab1VV_v(i+1) = update + 1;
  });
  end_gpu_timer(__KERNEL_NAME__);
  MBrr.dimensionner(n2,tab1RR(n2)-1);
  MBrv.dimensionner(n2,n1-n2,tab1RV(n2)-1);
  MBvv.dimensionner(n1-n2,n1-n2,tab1VV(n1-n2)-1);
 
  auto tab2RR_v = tab2RR.view_rw();
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(0, n2), KOKKOS_LAMBDA(const int i)
  {
    tab2RR_v(tab1RR_v(i) - 1) = i+1; // Diagonale
    rang_voisinRR_v(i) = tab1RR_v(i);
    rang_voisinRV_v(i) = tab1RV_v(i) - 1;
  });
  end_gpu_timer(__KERNEL_NAME__);
  auto tab2VV_v = tab2VV.view_rw();
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(0, n1-n2), KOKKOS_LAMBDA(const int i)
  {
    tab2VV_v(tab1VV_v(i) - 1) = i+1; // Diagonale
    rang_voisinVV_v(i) = tab1VV_v(i);
  });
  end_gpu_timer(__KERNEL_NAME__);
 
  int dim = Objet_U::dimension;
 
  MBrr.dimensionner(n2,tab1RR(n2)-1);
  MBrv.dimensionner(n2,n1-n2,tab1RV(n2)-1);
  MBvv.dimensionner(n1-n2,n1-n2,tab1VV(n1-n2)-1);
 
  auto tab2RV_v = tab2RV.view_rw();
  DoubleArrView coeffRR = MBrr.get_set_coeff().view_rw();
  DoubleArrView coeffRV = MBrv.get_set_coeff().view_rw();
  DoubleArrView coeffVV = MBvv.get_set_coeff().view_rw();
  CDoubleTabView face_normales = le_dom.face_normales().view_ro();
  CDoubleTabView inverse_quantitee_entrelacee_v = inverse_quantitee_entrelacee.view_ro();
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(ndeb, nfin), KOKKOS_LAMBDA(const int num_face)
  {
    int elem1 = face_voisins(num_face, 0);
    int elem2 = face_voisins(num_face, 1);
    const bool is_face_virt_bord = (num_face >= nb_faces) && (ind_faces_virt_bord(num_face - nb_faces) != -1);
    if (!is_face_virt_bord && elem1 != -1 && elem2 != -1)
      {
        double val = 0.;
        for (int d = 0; d < dim; d++)
          val += face_normales(num_face, d) * face_normales(num_face, d) * inverse_quantitee_entrelacee_v(num_face, d);
 
        // diagonale :
        if (elem1 < n2) Kokkos::atomic_add(&coeffRR(tab1RR_v(elem1) - 1), val);
        else            Kokkos::atomic_add(&coeffVV(tab1VV_v(elem1 - n2) - 1), val);
        if (elem2 < n2) Kokkos::atomic_add(&coeffRR(tab1RR_v(elem2) - 1), val);
        else            Kokkos::atomic_add(&coeffVV(tab1VV_v(elem2 - n2) - 1), val);
 
        if (elem1 > elem2)
          {
            if (elem1 < n2)
              {
                auto slot = Kokkos::atomic_fetch_add(&rang_voisinRR_v(elem2), 1);
                tab2RR_v(slot) = elem1 + 1;
                coeffRR(slot) -= val;
              }
            else
              {
                if (elem2 < n2)
                  {
                    auto slot = Kokkos::atomic_fetch_add(&rang_voisinRV_v(elem2), 1);
                    tab2RV_v(slot) = (elem1 - n2) + 1;
                    coeffRV(slot) -= val;
                  }
                else
                  {
                    auto slot = Kokkos::atomic_fetch_add(&rang_voisinVV_v(elem2 - n2), 1);
                    tab2VV_v(slot) = (elem1 - n2) + 1;
                    coeffVV(slot) -= val;
                  }
              }
          }
        else
          {
            if (elem2 < n2)
              {
                auto slot = Kokkos::atomic_fetch_add(&rang_voisinRR_v(elem1), 1);
                tab2RR_v(slot) = elem2 + 1;
                coeffRR(slot) -= val;
              }
            else
              {
                if (elem1 < n2)
                  {
                    auto slot = Kokkos::atomic_fetch_add(&rang_voisinRV_v(elem1), 1);
                    tab2RV_v(slot) = (elem2 - n2) + 1;
                    coeffRV(slot) -= val;
                  }
                else
                  {
                    auto slot = Kokkos::atomic_fetch_add(&rang_voisinVV_v(elem1 - n2), 1);
                    tab2VV_v(slot) = (elem2 - n2) + 1;
                    coeffVV(slot) -= val;
                  }
              }
          }
      }
  });
  end_gpu_timer(__KERNEL_NAME__);
  // On traite les conditions aux limites
  for (int i=0; i<les_cl.size(); i++)
    {
 
      // Le traitement depend du type de la condition aux limites :
      //  - Si condition de Neumann_sortie_libre ou Robin_VEF
      //  il faut calculer le coefficient sur la face et le prendre
      //  en compte dans le terme diagonal associe a l'element voisin.
      // - Si face de Cl avec une autre condition aux limites pas de
      // contribution a la matrice de pression.
 
      const Cond_lim& la_cl = les_cl[i];
      const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
      int nb_faces_bord_tot = le_bord.nb_faces_tot();
      if (sub_type(Neumann_sortie_libre,la_cl.valeur()) || sub_type(Robin_VEF,la_cl.valeur())  )
        {
          has_P_ref=1;
          MBrr.set_est_definie(1);
          CIntArrView front_num_face = le_bord.num_face().view_ro();
          DoubleArrView coeff_pression = static_cast<ArrOfDouble&>(les_coeff_pression).view_rw();
          const int coeff_pression_size = les_coeff_pression.size_array();
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(0, nb_faces_bord_tot), KOKKOS_LAMBDA(const int ind_face)
          {
            int num_face = front_num_face(ind_face);
            double val = 0.;
            for (int d = 0; d < dim; d++)
              val += face_normales(num_face, d) * face_normales(num_face, d) * inverse_quantitee_entrelacee_v(num_face, d);
 
            int elem = face_voisins(num_face, 0);
            if (elem < n2) Kokkos::atomic_add(&coeffRR(tab1RR_v(elem) - 1), val);
            else           Kokkos::atomic_add(&coeffVV(tab1VV_v(elem - n2) - 1), val);
            // On stocke les coefficients de pression sur les faces reelles
            if (num_face < coeff_pression_size)
              coeff_pression(num_face) = 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());
          CIntArrView front_num_face_coeff = le_bord.num_face().view_ro();
          CIntArrView face_associee_coeff = la_cl_perio.face_associee().view_ro();
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(0, nb_faces_bord_tot), KOKKOS_LAMBDA(const int ind_face)
          {
            if (ind_face < face_associee_coeff(ind_face)) // Process each periodic pair once
              {
                int num_face = front_num_face_coeff(ind_face);
                int elem1 = face_voisins(num_face, 0);
                int elem2 = face_voisins(num_face, 1);
                double val = 0.;
                for (int d = 0; d < dim; d++)
                  val += face_normales(num_face, d) * face_normales(num_face, d) * inverse_quantitee_entrelacee_v(num_face, d);
 
                // diagonale :
                if (elem1 < n2) Kokkos::atomic_add(&coeffRR(tab1RR_v(elem1) - 1), val);
                else            Kokkos::atomic_add(&coeffVV(tab1VV_v(elem1 - n2) - 1), val);
                if (elem2 < n2) Kokkos::atomic_add(&coeffRR(tab1RR_v(elem2) - 1), val);
                else            Kokkos::atomic_add(&coeffVV(tab1VV_v(elem2 - n2) - 1), val);
 
                if (elem1 > elem2)
                  {
                    if (elem1 < n2)
                      {
                        auto slot = Kokkos::atomic_fetch_add(&rang_voisinRR_v(elem2), 1);
                        tab2RR_v(slot) = elem1 + 1;
                        coeffRR(slot) -= val;
                      }
                    else
                      {
                        if (elem2 < n2)
                          {
                            auto slot = Kokkos::atomic_fetch_add(&rang_voisinRV_v(elem2), 1);
                            tab2RV_v(slot) = elem1 - n2 + 1;
                            coeffRV(slot) -= val;
                          }
                        else
                          {
                            auto slot = Kokkos::atomic_fetch_add(&rang_voisinVV_v(elem2 - n2), 1);
                            tab2VV_v(slot) = elem1 - n2 + 1;
                            coeffVV(slot) -= val;
                          }
                      }
                  }
                else
                  {
                    if (elem2 < n2)
                      {
                        auto slot = Kokkos::atomic_fetch_add(&rang_voisinRR_v(elem1), 1);
                        tab2RR_v(slot) = elem2 + 1;
                        coeffRR(slot) -= val;
                      }
                    else
                      {
                        if (elem1 < n2)
                          {
                            auto slot = Kokkos::atomic_fetch_add(&rang_voisinRV_v(elem1), 1);
                            tab2RV_v(slot) = elem2 - n2 + 1;
                            coeffRV(slot) -= val;
                          }
                        else
                          {
                            auto slot = Kokkos::atomic_fetch_add(&rang_voisinVV_v(elem1 - n2), 1);
                            tab2VV_v(slot) = elem2 - n2 + 1;
                            coeffVV(slot) -= val;
                          }
                      }
                  }
              }
          });
          end_gpu_timer(__KERNEL_NAME__);
        }
    }
  has_P_ref = (int)mp_max(has_P_ref);
  return 1;
}
 
/*! @brief Assemble la matrice de pression pour un fluide quasi compressible laplacein(P) est remplace par div(grad(P)/rho).
 *
 * @param (DoubleTab& tab_rho) mass volumique
 * @return (int) renvoie toujours 1
 */
int Assembleur_P_VEF::assembler_QC(const DoubleTab& tab_rho, Matrice& matrice)
{
  Cerr << "Assemblage de la matrice de pression pour Quasi Compressible en cours..." << finl;
  int stat=assembler(matrice);
  set_resoudre_en_u(0);
  return stat;
}
 
int Assembleur_P_VEF::modifier_secmem(DoubleTab& secmem)
{
  const Domaine_VEF& le_dom = le_dom_VEF.valeur();
  const Domaine_Cl_VEF& le_dom_cl = le_dom_Cl_VEF.valeur();
  int nb_cond_lim = le_dom_cl.nb_cond_lim();
  const IntTab& face_voisins = le_dom.face_voisins();
 
  // Modification du second membre :
  for (int i=0; i<nb_cond_lim; i++)
    {
      const Cond_lim_base& la_cl_base = le_dom_cl.les_conditions_limites(i).valeur();
      const Front_VF& la_front_dis = ref_cast(Front_VF,la_cl_base.frontiere_dis());
      const Champ_front_base& champ_front = la_cl_base.champ_front();
      int ndeb = la_front_dis.num_premiere_face();
      int nfin = ndeb + la_front_dis.nb_faces();
 
      // GF on est passe en increment de pression
      if ((sub_type(Neumann_sortie_libre,la_cl_base)) && (!get_resoudre_increment_pression()))
        {
          const Neumann_sortie_libre& la_cl_Neumann = ref_cast(Neumann_sortie_libre, la_cl_base);
          ToDo_Kokkos("critical");
          for (int num_face=ndeb; num_face<nfin; num_face++)
            {
              double Pimp = la_cl_Neumann.flux_impose(num_face-ndeb);
              double coef = les_coeff_pression[num_face]*Pimp;
              secmem[face_voisins(num_face,0)] += coef;
            }
        }
 
      /*if ((sub_type(Robin_VEF,la_cl_base)) )
        {
          // [vkr/oswr] we should modify the value of increment_pression_bord in order to respect the oswr algorithm
          double Pstar_OSWR, coef;
          const Robin_VEF& la_cl_robin = ref_cast(Robin_VEF, la_cl_base);
          for (int num_face=ndeb; num_face<nfin; num_face++)
            {
              Pstar_OSWR = la_cl_robin.increment_pression_bord(num_face);
              coef = les_coeff_pression[num_face]*Pstar_OSWR+100;
              Cerr << "PASSAGE PAR GET RESOUDRE INCR PRESSION" <<finl;
              secmem[face_voisins(num_face,0)] += coef;
            }
        }*/
 
      else if ( champ_front.instationnaire() && get_resoudre_en_u() )
        {
          const DoubleTab& Gpt = champ_front.derivee_en_temps();
          bool ch_unif = (Gpt.nb_dim()==1);
          ToDo_Kokkos("critical");
          for (int num_face=ndeb; num_face<nfin; num_face++)
            {
              double Stt = 0.;
              for (int k=0; k<dimension; k++)
                {
                  double Gpoint = ch_unif ? Gpt(k) : Gpt(num_face - ndeb, k);
                  Stt -= Gpoint * le_dom.face_normales(num_face, k);
                }
              secmem(face_voisins(num_face,0)) += Stt;
            }
        }
    }
  secmem.echange_espace_virtuel();
  return 1;
}
 
int Assembleur_P_VEF::modifier_solution(DoubleTab& pression)
{
  // Projection :
  double press_0;
  if(!has_P_ref)
    {
      // On prend la pression minimale comme pression de reference
      // afin d'avoir la meme pression de reference en sequentiel et parallele
      press_0=DMAXFLOAT;
      int nb_elem=le_dom_VEF->domaine().nb_elem();
      ToDo_Kokkos("critical");
      for(int n=0; n<nb_elem; n++)
        if (pression[n] < press_0)
          press_0 = pression[n];
      press_0 = Process::mp_min(press_0);
      ToDo_Kokkos("critical");
      for(int n=0; n<nb_elem; n++)
        pression[n] -=press_0;
 
      pression.echange_espace_virtuel();
    }
  return 1;
}
 
/*! @brief Modifier eventuellement la matrice pour la rendre definie si elle ne l'est pas Valeurs par defaut:
 *
 *     Contraintes:
 *     Acces: entree
 *
 * @return (int) renvoie 1 si la matrice est modifiee 0 sinon
 */
int Assembleur_P_VEF::modifier_matrice(Matrice& matrice)
{
  int matrice_modifiee=0;
  Matrice_Bloc& mat_bloc = ref_cast(Matrice_Bloc, matrice.valeur());
  Matrice_Morse_Sym& A00RR = ref_cast(Matrice_Morse_Sym,mat_bloc.get_bloc(0,0).valeur());
  // Recherche de l'element sur lequel on impose la pression de reference
  const bool is_first_proc_with_real_elems = Process::me() == Process::mp_min(le_dom_VEF->nb_elem() ? Process::me() : 1e8);
  if (is_first_proc_with_real_elems && !A00RR.get_est_definie())
    {
      int element_referent=0;
      double distance=DMAXFLOAT;
      const DoubleTab& coord=le_dom_VEF->xp();
      int n = le_dom_VEF->nb_elem();
      ToDo_Kokkos("critical");
      for(int i=0; i<n; i++)
        {
          double tmp=0;
          for (int j=0; j<dimension; j++)
            tmp+=coord(i,j)*coord(i,j);
          if (inf_strict(tmp,distance) && !est_egal(A00RR(i,i),0.))
            {
              distance=tmp;
              element_referent=i;
            }
        }
      Cerr << "On modifie la ligne (element) " << element_referent << finl;
      A00RR(element_referent,element_referent)*=2;
      matrice_modifiee=1;
    }
  // has_P_ref=1;
  A00RR.set_est_definie(1);
  return matrice_modifiee;
}
 
const Domaine_dis_base& Assembleur_P_VEF::domaine_dis_base() const
{
  return le_dom_VEF.valeur();
}
 
const Domaine_Cl_dis_base& Assembleur_P_VEF::domaine_Cl_dis_base() const
{
  return le_dom_Cl_VEF.valeur();
}
 
void Assembleur_P_VEF::associer_domaine_dis_base(const Domaine_dis_base& le_dom_dis)
{
  le_dom_VEF = ref_cast(Domaine_VEF,le_dom_dis);
}
 
void Assembleur_P_VEF::associer_domaine_cl_dis_base(const Domaine_Cl_dis_base& le_dom_Cl_dis)
{
  le_dom_Cl_VEF = ref_cast(Domaine_Cl_VEF, le_dom_Cl_dis);
}
 
void Assembleur_P_VEF::completer(const Equation_base& Eqn)
{
  mon_equation=Eqn;
}