Assembleur_P_VEFPreP1B.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_VEFPreP1B.h>
#include <Matrice_Bloc_Sym.h>
#include <Domaine.h>
#include <Dirichlet.h>
#include <Dirichlet_homogene.h>
#include <Dirichlet_entree_fluide_leaves.h>
#include <Neumann_sortie_libre.h>
#include <Milieu_base.h>
#include <Navier_Stokes_std.h>
#include <TRUSTLists.h>
#include <Op_Grad_VEF_P1B_Face.h>
#include <Op_Div_VEFP1B_Elem.h>
#include <SFichier.h>
#include <Champ_Fonc_P1NC.h>
#include <SolveurPP1B.h>
#include <Check_espace_virtuel.h>
#include <Solv_Petsc.h>
#include <Solv_AMG.h>
#include <Matrice_Petsc.h>
#include <Robin_VEF.h>
 
Implemente_instanciable(Assembleur_P_VEFPreP1B,"Assembleur_P_VEFPreP1B",Assembleur_P_VEF);
 
// printOn and readOn
 
Sortie& Assembleur_P_VEFPreP1B::printOn(Sortie& s ) const
{
  return Assembleur_base::printOn(s);
}
 
Entree& Assembleur_P_VEFPreP1B::readOn(Entree& is )
{
  return Assembleur_base::readOn(is);
}
 
static double estim_cond(Matrice& A)
{
  int i=A.valeur( ).nb_lignes();
  //Cout << "rang " << i << finl;
  DoubleVect p(A.valeur( ).nb_lignes());
  for(int j=0; j<i; j++)
    p(j)=drand48();
  double normep=sqrt(mp_prodscal(p,p));
  p/=normep;
  DoubleVect r(A.valeur( ).nb_lignes());
  double lmax=1.e-36, lmin=1.e24;
  //double cond=-1;
  i=10000;
  while(i-- )
    {
      //cond=lmax/lmin;
      A.valeur( ).multvect(p,r);
      double pscpr=mp_prodscal(p,r);
      //Cout << "lambda " << i << " : " << pscpr << " , " << lmax/lmin << finl;
      double pscrr=mp_prodscal(r,r);
      double epsilon=-pscpr/pscrr;
      p.ajoute(epsilon,r);
      double normepbis=sqrt(mp_prodscal(p,p));
      //Cout << "normep : " << normep << finl;
      p/=normepbis;
      if(pscpr>lmax)
        lmax=pscpr;
      else if(pscpr<lmin)
        lmin=pscpr;
    }
  //Cout << "lmax " << lmax << " lmin = " << lmin << finl;
  //Cout << p << finl;
  if (A.valeur( ).nb_lignes()<30) A.valeur( ).imprimer_formatte(Cerr);
 
  return(lmax/lmin);
}
 
void Assembleur_P_VEFPreP1B::completer(const Equation_base& eqn)
{
  mon_equation=eqn;
  Assembleur_P_VEF::completer(eqn);
  if (domaine_Vef().get_P1Bulle())
    {
      // To change basis and recover P1Bulle
      alpha_=1./Objet_U::dimension;
      beta_=1./(Objet_U::dimension*(Objet_U::dimension+1));
    }
}
 
const Domaine_VEF& Assembleur_P_VEFPreP1B::domaine_Vef() const
{
  return ref_cast(Domaine_VEF, le_dom_VEF.valeur());
}
 
extern void assemblerP0P0(const Domaine_dis_base& z,
                          const Domaine_Cl_dis_base& zcl,
                          Matrice& matrice,
                          const DoubleTab& inverse_quantitee_entrelacee);
 
extern void assemblerP1P1(const Domaine_dis_base& z,
                          const Domaine_Cl_dis_base& zcl,
                          Matrice& matrice,
                          const DoubleTab& inverse_quantitee_entrelacee, const ArrOfDouble& coef_som);
 
extern void assemblerPaPa(const Domaine_dis_base& z,
                          const Domaine_Cl_dis_base& zcl,
                          Matrice& matrice, const DoubleTab& inverse_quantitee_entrelacee);
 
extern void assemblerP0P1(const Domaine_dis_base& z,
                          const Domaine_Cl_dis_base& zcl,
                          Matrice& matrice, const DoubleTab& inverse_quantitee_entrelacee, const ArrOfDouble& coef_som);
 
extern void assemblerP0Pa(const Domaine_dis_base& z,
                          const Domaine_Cl_dis_base& zcl,
                          Matrice& matrice,
                          const DoubleTab& inverse_quantitee_entrelacee);
 
extern  void assemblerP1Pa(const Domaine_dis_base& z,
                           const Domaine_Cl_dis_base& zcl,
                           Matrice& matrice,
                           const DoubleTab& inverse_quantitee_entrelacee, const ArrOfDouble& coef_som);
 
extern void updateP0P0(const Domaine_dis_base& z,
                       const Domaine_Cl_dis_base& zcl,
                       Matrice& matrice,
                       const DoubleTab& inverse_quantitee_entrelacee);
 
extern void updateP1P1(const Domaine_dis_base& z,
                       const Domaine_Cl_dis_base& zcl,
                       Matrice& matrice,
                       const DoubleTab& inverse_quantitee_entrelacee, const ArrOfDouble& coef_som);
extern void modifieP1P1neumann(const Domaine_dis_base& z,
                               const Domaine_Cl_dis_base& zcl,
                               Matrice& matrice,
                               const DoubleTab& inverse_quantitee_entrelacee, const ArrOfDouble& coef_som);
extern void updatePaPa(const Domaine_dis_base& z,
                       const Domaine_Cl_dis_base& zcl,
                       Matrice& matrice,
                       const DoubleTab& inverse_quantitee_entrelacee);
 
extern void updateP0P1(const Domaine_dis_base& z,
                       const Domaine_Cl_dis_base& zcl,
                       Matrice& matrice,
                       const DoubleTab& inverse_quantitee_entrelacee, const ArrOfDouble& coef_som);
 
extern void updateP0Pa(const Domaine_dis_base& z,
                       const Domaine_Cl_dis_base& zcl,
                       Matrice& matrice,
                       const DoubleTab& inverse_quantitee_entrelacee);
 
extern  void updateP1Pa(const Domaine_dis_base& z,
                        const Domaine_Cl_dis_base& zcl,
                        Matrice& matrice,
                        const DoubleTab& inverse_quantitee_entrelacee, const ArrOfDouble& coef_som);
 
extern int verifier( const Assembleur_P_VEFPreP1B& ass,
                     const Matrice_Bloc_Sym& matrice,
                     const Domaine_VEF& domaine_VEF,
                     const DoubleTab& inverse_quantitee_entrelacee);
 
int Assembleur_P_VEFPreP1B::assembler(Matrice& la_matrice)
{
  return Assembleur_P_VEF::assembler(la_matrice);
}
 
int Assembleur_P_VEFPreP1B::assembler_rho_variable(Matrice& la_matrice, const Champ_Don_base& rho)
{
  // Multiply by mass density at faces
  if (!sub_type(Champ_Fonc_P1NC, rho))
    {
      Cerr << "The mass density is not at faces in Assembleur_P_VEFPreP1B::assembler_rho_variable." << finl;
      Process::exit();
    }
  const DoubleVect& volumes_entrelaces=domaine_Vef().volumes_entrelaces();
  const DoubleVect& masse_volumique=rho.valeurs();
  DoubleVect quantitee_entrelacee(volumes_entrelaces);
  int size=quantitee_entrelacee.size_array();
  for (int i=0; i<size; i++)
    quantitee_entrelacee(i)=(volumes_entrelaces(i)*masse_volumique(i));
 
  // Assemble the matrix
  return assembler_mat(la_matrice,quantitee_entrelacee,1,1);
}
 
void zero(Matrice_Bloc_Sym& matrice)
{
  int nombre_supports = matrice.ordre();
  for(int i=0; i<nombre_supports; i++)
    for(int j=i; j<nombre_supports; j++)
      {
        Matrice_Bloc& bloc_ij=ref_cast(Matrice_Bloc, matrice.get_bloc(i,j).valeur());
        ref_cast(Matrice_Morse, bloc_ij.get_bloc(0,0).valeur()).clean();
        if(Process::is_parallel())
          {
            ref_cast(Matrice_Morse, bloc_ij.get_bloc(0,1).valeur()).clean();
            ref_cast(Matrice_Morse, bloc_ij.get_bloc(1,0).valeur()).clean();
            ref_cast(Matrice_Morse, bloc_ij.get_bloc(1,1).valeur()).clean();
          }
      }
}
 
 
 
int Assembleur_P_VEFPreP1B::assembler_mat(Matrice& la_matrice,const DoubleVect& quantitee_entrelacee, int incr_pression, int resoudre_en_u)
{
  // Set the flags of Assembleur_base
  set_resoudre_increment_pression(incr_pression);
  set_resoudre_en_u(resoudre_en_u);
 
  SolveurSys& solveur_pression = ref_cast(Navier_Stokes_std, mon_equation.valeur()).solveur_pression();
  if (solveur_pression->read_matrix()) // Read the matrix from a file (only supported by PETSc for now)
    {
      la_matrice.typer("Matrice_Petsc");
    }
  else // Assemble the matrix
    {
      const Domaine_VEF& domaine_vef = domaine_Vef();
      Cerr << "Assembling pressure matrix" << (domaine_vef.get_alphaE() ? " P0" : "")
           << (domaine_vef.get_alphaS() ? " P1" : "") << (domaine_vef.get_alphaA() ? " Pa" : "") << " in progress..." << finl;
 
      // Partitions must have ghost cell width of 2
      // if P1 or Pa support is used...
      if (Process::is_parallel() &&
          domaine_vef.domaine().nb_joints() &&
          domaine_vef.domaine().joint(0).epaisseur() < 2 &&
          (domaine_vef.get_alphaS() || domaine_vef.get_alphaA()))
        {
          Cerr << "Ghost cells width of " << domaine_vef.domaine().joint(0).epaisseur()
               << " is not enough for assembling VEFPreP1B pressure matrix" << finl;
          Cerr << "for parallel calculation. Partition your mesh with larg_joint option set to 2." << finl;
          Process::exit();
        }
 
      DoubleTab inverse_quantitee_entrelacee;
      const Domaine_Cl_VEF& domaine_Cl_VEF = le_dom_Cl_VEF.valeur();
      calculer_inv_volume(inverse_quantitee_entrelacee, domaine_Cl_VEF, quantitee_entrelacee);
      int P0 = 0;
      int P1 = P0 + domaine_vef.get_alphaE();
      int Pa = P1 + domaine_vef.get_alphaS();
      int nombre_supports = Pa + domaine_vef.get_alphaA();
      DoubleVect coef_som;
      if (domaine_vef.get_alphaS())
        {
          domaine_vef.domaine().creer_tableau_elements(coef_som);
          for (int elem = 0; elem < coef_som.size_totale(); elem++)
            coef_som[elem] = Op_Grad_VEF_P1B_Face::calculer_coef_som(elem, domaine_Cl_VEF, domaine_vef);
          coef_som.echange_espace_virtuel();
        }
 
      // Assemble the complete matrix according to the chosen supports
      Matrice la_matrice_de_travail_;                // Working matrix
      la_matrice_de_travail_.typer("Matrice_Bloc_Sym");
      Matrice_Bloc_Sym& la_matrice_bloc_sym_de_travail = ref_cast(Matrice_Bloc_Sym, la_matrice_de_travail_.valeur());
      if (la_matrice_bloc_sym_de_travail.nb_bloc_lignes()==0)
        {
          la_matrice_bloc_sym_de_travail.dimensionner(nombre_supports, nombre_supports);
          if (domaine_vef.get_alphaE())
            assemblerP0P0(domaine_vef, domaine_Cl_VEF, la_matrice_bloc_sym_de_travail.get_bloc(P0, P0),
                          inverse_quantitee_entrelacee);
 
          if (domaine_vef.get_alphaS())
            assemblerP1P1(domaine_vef, domaine_Cl_VEF, la_matrice_bloc_sym_de_travail.get_bloc(P1, P1),
                          inverse_quantitee_entrelacee, coef_som);
 
          if (domaine_vef.get_alphaE() && domaine_vef.get_alphaS())
            assemblerP0P1(domaine_vef, domaine_Cl_VEF, la_matrice_bloc_sym_de_travail.get_bloc(P0, P1),
                          inverse_quantitee_entrelacee, coef_som);
 
          if (domaine_vef.get_alphaA())
            assemblerPaPa(domaine_vef, domaine_Cl_VEF, la_matrice_bloc_sym_de_travail.get_bloc(Pa, Pa),
                          inverse_quantitee_entrelacee);
 
          if (domaine_vef.get_alphaS() && domaine_vef.get_alphaA())
            assemblerP1Pa(domaine_vef, domaine_Cl_VEF, la_matrice_bloc_sym_de_travail.get_bloc(P1, Pa),
                          inverse_quantitee_entrelacee, coef_som);
 
          if (domaine_vef.get_alphaE() && domaine_vef.get_alphaA())
            assemblerP0Pa(domaine_vef, domaine_Cl_VEF, la_matrice_bloc_sym_de_travail.get_bloc(P0, Pa),
                          inverse_quantitee_entrelacee);
        }
      // Set the matrix to zero even if it was correctly built and filled by the
      // assemblerPiPi methods, then refill it. Why? To have sufficient test coverage
      // of the updatePiPi methods while waiting to properly factorize the
      // assemblerPiPi and updatePiPi methods
      zero(la_matrice_bloc_sym_de_travail);
      {
        if (domaine_vef.get_alphaE())
          updateP0P0(domaine_vef, domaine_Cl_VEF, la_matrice_bloc_sym_de_travail.get_bloc(P0, P0),
                     inverse_quantitee_entrelacee);
 
        if (domaine_vef.get_alphaS())
          {
            updateP1P1(domaine_vef, domaine_Cl_VEF, la_matrice_bloc_sym_de_travail.get_bloc(P1, P1),
                       inverse_quantitee_entrelacee, coef_som);
            if (domaine_vef.get_cl_pression_sommet_faible() == 0)
              modifieP1P1neumann(domaine_vef, domaine_Cl_VEF, la_matrice_bloc_sym_de_travail.get_bloc(P1, P1),
                                 inverse_quantitee_entrelacee, coef_som);
          }
        if (domaine_vef.get_alphaE() && domaine_vef.get_alphaS())
          updateP0P1(domaine_vef, domaine_Cl_VEF, la_matrice_bloc_sym_de_travail.get_bloc(P0, P1),
                     inverse_quantitee_entrelacee, coef_som);
 
        if (domaine_vef.get_alphaA())
          updatePaPa(domaine_vef, domaine_Cl_VEF, la_matrice_bloc_sym_de_travail.get_bloc(Pa, Pa),
                     inverse_quantitee_entrelacee);
 
        if (domaine_vef.get_alphaS() && domaine_vef.get_alphaA())
          updateP1Pa(domaine_vef, domaine_Cl_VEF, la_matrice_bloc_sym_de_travail.get_bloc(P1, Pa),
                     inverse_quantitee_entrelacee, coef_som);
 
        if (domaine_vef.get_alphaE() && domaine_vef.get_alphaA())
          updateP0Pa(domaine_vef, domaine_Cl_VEF, la_matrice_bloc_sym_de_travail.get_bloc(P0, Pa),
                     inverse_quantitee_entrelacee);
      }
      //trustIdType ordre_matrice = mp_sum(la_matrice_bloc_sym_de_travail.nb_lignes());
      //Cerr << "Order of the matrix = " << ordre_matrice << finl;
 
      // Verification method
      char *theValue = getenv("TRUST_VERIFIE_MATRICE_VEF");
      if (theValue != nullptr) verifier(*this, la_matrice_bloc_sym_de_travail, domaine_vef, inverse_quantitee_entrelacee);
 
      ////////////////////////////////////////////
      // Optional basis change P0P1->P1Bulle
      ////////////////////////////////////////////
      if (changement_base())
        changer_base_matrice(la_matrice_de_travail_); // A->A~
 
      /////////////////////////////////////////////////////////////
      // Modify the matrix if there is no reference pressure
      /////////////////////////////////////////////////////////////
      modifier_matrice(la_matrice_de_travail_);
 
      // Optional conversion to Matrice_Morse_Sym
      if (ref_cast(Navier_Stokes_std, mon_equation.valeur()).solveur_pression()->supporte_matrice_morse_sym() &&
          domaine_vef.get_alphaE() != nombre_supports) // We are not in P0 only
        {
          //////////////////////////////////////////////////////////
          // The returned matrix is a cleaned Matrice_Morse_Sym
          // if the solver used supports this matrix type
          // and if we are not in P0 only (in that case,
          // conversion to Mat_Morse_Sym brings nothing and crashes SSOR
          // because the Morse matrix then contains virtual parts
          // with no shared items and we hit
          // the assert assert(*tab2_ptr<=n); Two other solutions:
          // -Modify SSOR to solve only the real part
          // -For P0 only, empty VV and RV in Mat_Bloc_Sym
          //////////////////////////////////////////////////////////
          la_matrice.typer("Matrice_Morse_Sym");
          ref_cast(Matrice_Bloc_Sym, la_matrice_de_travail_.valeur()).BlocSymToMatMorseSym(
            ref_cast(Matrice_Morse_Sym, la_matrice.valeur()));
          ref_cast(Matrice_Morse_Sym, la_matrice.valeur()).compacte(
            2);// Removal of zero and near-zero coefficients
        }
      else
        {
          /////////////////////////////////////////////////////////
          // The returned matrix is a cleaned Matrice_Bloc_Sym
          /////////////////////////////////////////////////////////
          la_matrice = la_matrice_de_travail_;
          for (int i = 0; i < nombre_supports; i++)
            for (int j = i; j < nombre_supports; j++)
              {
                Matrice_Bloc_Sym& mat_bloc_sym = ref_cast(Matrice_Bloc_Sym, la_matrice.valeur());
                Matrice_Bloc& bloc_ij = ref_cast(Matrice_Bloc, mat_bloc_sym.get_bloc(i, j).valeur());
                ref_cast(Matrice_Morse, bloc_ij.get_bloc(0, 0).valeur()).compacte(
                  2); // Removal of zero and near-zero coefficients
              }
        }
    }
 
  // If not already done, use a solver (SolveurPP1B) that will perform basis changes for the solution and the right-hand side
  if (changement_base() && solveur_pression->que_suis_je() != "SolveurPP1B")
    {
      SolveurSys solveur_pression_lu = solveur_pression;
      solveur_pression.typer("SolveurPP1B");
      SolveurPP1B& solveur_pression_PP1B = ref_cast(SolveurPP1B, solveur_pression.valeur());
      solveur_pression_PP1B.associer(*this, solveur_pression_lu);
    }
 
  //////////////////////////////////////////////////////
  // Optional display of the matrix condition number
  //////////////////////////////////////////////////////
  char* theValue2 = getenv("TRUST_CONDITIONNEMENT_MATRICE");
  if(theValue2 != nullptr)
    Cout << "Estimation of the matrix condition number: " << estim_cond(la_matrice)<<finl;
 
  return 1;
}
 
int Assembleur_P_VEFPreP1B::modifier_secmem(DoubleTab& b)
{
  const Domaine_VEF& le_dom = domaine_Vef();
 
  // Verification on the Pa support
  if (le_dom.get_alphaA())
    {
      // Verify on edges that:
      const ArrOfInt& renum_arete_perio=le_dom.get_renum_arete_perio();
      const ArrOfInt& ok_arete=le_dom.get_ok_arete();
      int npa=le_dom.numero_premiere_arete();
      // b does not necessarily have its virtual space up to date
      int nb_aretes=le_dom.domaine().nb_aretes();
      ToDo_Kokkos("critical");
      for(int i=0; i<nb_aretes; i++)
        if(!ok_arete[i] && b(npa+i)!=0.) // Superfluous edges have a zero value
          {
            Cerr << "Pb div Edges, the pressure on edge " << i << " (which is superfluous) is not zero." << finl;
            Process::exit();
          }
        else if ( (renum_arete_perio[i]!=i) && b(npa+i)!=0.) // Periodic edges have a zero value
          {
            Cerr << "Pb div Periodic Edges, the pressure on edge " << i << " (which is periodic) is not zero." << finl;
            Process::exit();
          }
    }
 
  // Verification on the P1 support
  if (le_dom.get_alphaS())
    {
      // Verify that the pressure on periodic vertices is zero
      const Domaine& dom=le_dom.domaine();
      int nps=le_dom.numero_premier_sommet();
      int ns=le_dom.domaine().nb_som();
      CIntArrView renum_som_perio = dom.get_renum_som_perio().view_ro();
      CDoubleArrView b_v = static_cast<const DoubleVect&>(b).view_ro();
      Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), ns, KOKKOS_LAMBDA(
                             const int i)
      {
        int k=renum_som_perio(i);
        if((k!=i)&& b_v(nps+i)!=0.)
          {
            printf("Pb div Som, the pressure on vertex %ld (which is periodic) is not zero.\n",(long)i);
            printf("In plain terms, the right-hand side is not zero on a periodic vertex.\n");
            if (b_v(nps+i)!=b_v(nps+k))
              {
                printf("Moreover, the right-hand side does not have the same value on the 2 periodic vertices.\n");
                printf("b(nps+i)=%f <> b(nps+k)=%f\n",b_v(nps+i),b_v(nps+k));
              }
            Kokkos::abort("There is a probability that the model used is poorly implemented for\na periodicity condition. Contact TRUST support.");
          }
      });
      end_gpu_timer(__KERNEL_NAME__);
    }
 
  if (get_resoudre_en_u())
    {
      const Domaine_VEF& domaine_VEF =  domaine_Vef();
      const Domaine_Cl_VEF& domaine_Cl = le_dom_Cl_VEF.valeur();
 
      CDoubleArrView porosite_face = domaine_Cl.equation().milieu().porosite_face().view_ro();
 
      const int nb_cond_lim = domaine_Cl.nb_cond_lim();
 
      /**************************/
      /* Retrieval of Gpoint */
      /**************************/
      DoubleTrav tab_Gpoint(equation().inconnue().valeurs());
      DoubleTabView Gpoint = tab_Gpoint.view_wo();
 
      //Gpoint=0.; A DoubleTrav is initialized to 0
      int Gpoint_nul = 1; // Flag to potentially save a echange_espace_virtuel call
      for (int cond_lim=0; cond_lim<nb_cond_lim; cond_lim++)
        {
          const Cond_lim_base& cl_base = domaine_Cl.les_conditions_limites(cond_lim).valeur();
 
          const Front_VF& front_VF = ref_cast(Front_VF,cl_base.frontiere_dis());
          const Champ_front_base& champ_front = cl_base.champ_front();
 
          const int ndeb = front_VF.num_premiere_face();
          const int nfin = ndeb+front_VF.nb_faces();
 
          /* Test on the nature of the field at the domain boundary */
          if (sub_type(Entree_fluide_vitesse_imposee, cl_base) && champ_front.instationnaire())
            {
              Gpoint_nul = 0;
              CDoubleTabView gpoint = champ_front.derivee_en_temps().view_ro();
              CDoubleArrView arr_gpoint = static_cast<const DoubleVect&>(champ_front.derivee_en_temps()).view_ro();
              bool ch_unif = (champ_front.derivee_en_temps().nb_dim()==1);
              int dimension_ = Objet_U::dimension;
              Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), Kokkos::RangePolicy<>(ndeb, nfin),
                                   KOKKOS_LAMBDA(const int num_face)
              {
                for (int dim=0; dim<dimension_; dim++)
                  Gpoint(num_face,dim) = porosite_face(num_face) * (ch_unif ? arr_gpoint(dim) : gpoint(num_face-ndeb,dim));
              });
              end_gpu_timer(__KERNEL_NAME__);
            }
        }
 
      //For parallelism
      if (!Gpoint_nul) tab_Gpoint.echange_espace_virtuel();
 
      /******************************/
      /* End of Gpoint retrieval */
      /******************************/
 
      /*********************************/
      /* Modification of the right-hand side */
      /*********************************/
 
      if (domaine_VEF.get_alphaE()) modifier_secmem_elem(tab_Gpoint,b);
      if (domaine_VEF.get_alphaS()) modifier_secmem_som(tab_Gpoint,b);
      if (domaine_VEF.get_alphaA()) modifier_secmem_aretes(tab_Gpoint,b);
 
      /**********************************/
      /* End of right-hand side modification */
      /**********************************/
 
      b.echange_espace_virtuel();
    }
 
  return 1;
}
 
int Assembleur_P_VEFPreP1B::modifier_secmem_elem(const DoubleTab& tab_Gpoint, DoubleTab& tab_b)
{
  const Domaine_VEF& domaine_VEF =  domaine_Vef();
  const Domaine_Cl_VEF& domaine_Cl = le_dom_Cl_VEF.valeur();
  const int nb_cond_lim = domaine_Cl.nb_cond_lim();
  for (int cond_lim=0; cond_lim<nb_cond_lim; cond_lim++)
    {
      const Cond_lim_base& cl_base = domaine_Cl.les_conditions_limites(cond_lim).valeur();
 
      const Front_VF& front_VF = ref_cast(Front_VF,cl_base.frontiere_dis());
      const Champ_front_base& champ_front = cl_base.champ_front();
      /* Test on the nature of the field at the domain boundary */
      if (sub_type(Entree_fluide_vitesse_imposee, cl_base) && champ_front.instationnaire() )
        {
          // Build the list of faces to process (real + virtual)
          const int nb_faces_bord_tot = front_VF.nb_faces_tot();
          int dimension_ = Objet_U::dimension;
          CIntTabView face_voisins = domaine_VEF.face_voisins().view_ro();
          CDoubleTabView face_normales = domaine_VEF.face_normales().view_ro();
          CDoubleTabView Gpoint = tab_Gpoint.view_ro();
          CIntArrView front_num_face = front_VF.num_face().view_ro();
          DoubleArrView b = static_cast<DoubleVect&>(tab_b).view_rw();
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_faces_bord_tot, KOKKOS_LAMBDA(const int ind_face)
          {
            const int num_face = front_num_face(ind_face);
            const int elem = face_voisins(num_face,0);
            assert(elem!=-1);
 
            for (int dim=0; dim<dimension_; dim++)
              Kokkos::atomic_add(&b(elem), -Gpoint(num_face,dim)*face_normales(num_face,dim));
          });
          end_gpu_timer(__KERNEL_NAME__);
        }
    }
 
  return 1;
}
 
int Assembleur_P_VEFPreP1B::modifier_secmem_som(const DoubleTab& tab_Gpoint, DoubleTab& tab_b)
{
  const Domaine_VEF& domaine_VEF =  domaine_Vef();
  const Domaine_Cl_VEF& domaine_Cl = le_dom_Cl_VEF.valeur();
  const Domaine& domaine = domaine_VEF.domaine();
 
  const int nb_cond_lim = domaine_Cl.nb_cond_lim();
  const int nb_elem_tot = (domaine_VEF.get_alphaE()?domaine.nb_elem_tot():0);
  const int nb_faces_elem = domaine.nb_faces_elem();
 
  const double coeff_dim = Objet_U::dimension*(Objet_U::dimension+1);
 
  for (int cond_lim=0; cond_lim<nb_cond_lim; cond_lim++)
    {
      const Cond_lim_base& cl_base = domaine_Cl.les_conditions_limites(cond_lim).valeur();
 
      const Front_VF& front_VF = ref_cast(Front_VF,cl_base.frontiere_dis());
      const Champ_front_base& champ_front = cl_base.champ_front();
 
      /* Test on the nature of the field at the domain boundary */
      if (sub_type(Entree_fluide_vitesse_imposee, cl_base)  && champ_front.instationnaire())
        {
          // Build the list of faces to process (real + virtual)
          const int nb_faces_bord_tot = front_VF.nb_faces_tot();
          int dimension_ = Objet_U::dimension;
          CIntTabView elem_faces = domaine_VEF.elem_faces().view_ro();
          CIntTabView face_sommets = domaine_VEF.face_sommets().view_ro();
          CIntTabView face_voisins = domaine_VEF.face_voisins().view_ro();
          CIntTabView elem_sommets = domaine.les_elems().view_ro();
          CDoubleTabView face_normales = domaine_VEF.face_normales().view_ro();
          CDoubleTabView Gpoint = tab_Gpoint.view_ro();
          CIntArrView renum_som_perio = domaine.get_renum_som_perio().view_ro();
          CIntArrView front_num_face = front_VF.num_face().view_ro();
          DoubleArrView b = static_cast<DoubleVect&>(tab_b).view_rw();
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_faces_bord_tot, KOKKOS_LAMBDA(const int ind_face)
          {
            double sigma[3];
            const int num_face = front_num_face(ind_face);
            const int elem = face_voisins(num_face,0);
            assert(elem!=-1);
 
            //Compute the velocity at the element centre
            for (int i_sigma = 0; i_sigma < 3; ++i_sigma)
              sigma[i_sigma] = 0;
 
            for (int face_loc=0; face_loc<nb_faces_elem; face_loc++)
              {
                const int face = elem_faces(elem,face_loc);
 
                for(int comp=0; comp<dimension_; comp++)
                  sigma[comp] += Gpoint(face, comp);
              }
 
            //Compute the divergence of the velocity
            for(int face_loc=0; face_loc<nb_faces_elem; face_loc++)
              {
                const int som = nb_elem_tot + renum_som_perio(elem_sommets(elem,face_loc));
                const int face = elem_faces(elem,face_loc);
 
                double psc=0;
                double signe=1.;
                if(elem != face_voisins(face,0)) signe=-1.;
 
                for(int comp=0; comp<dimension_; comp++)
                  psc += sigma[comp]*face_normales(face,comp);
 
                Kokkos::atomic_add(&b(som), -signe*psc/coeff_dim);
              }
 
            double flux = 0. ;
            for (int comp=0; comp<dimension_; comp++)
              flux += Gpoint(num_face,comp) * face_normales(num_face,comp) ;
 
            flux*=1./dimension_;
            for(int som_loc=0; som_loc<nb_faces_elem-1; som_loc++)
              {
                const int som = nb_elem_tot + renum_som_perio(face_sommets(num_face,som_loc));
                // DOUBT_HARI: Should it be atomic?
                Kokkos::atomic_add(&b(som), -flux);
              }
            //End of velocity divergence computation
          });
          end_gpu_timer(__KERNEL_NAME__);
        }
    }
 
  return 1;
}
 
int Assembleur_P_VEFPreP1B::modifier_secmem_aretes(const DoubleTab& Gpoint, DoubleTab& b)
{
  return 1;
}
 
int Assembleur_P_VEFPreP1B::modifier_solution(DoubleTab& tab_pression)
{
 
  //  if (!has_P_ref) exit();
  const Domaine_VEF& le_dom = domaine_Vef();
 
  // Verification on edges
  if (le_dom.get_alphaA())
    {
      // Impose pressure = 0 on superfluous edges:
      const IntVect& ok_arete=le_dom.get_ok_arete();
      const ArrOfInt& renum_arete_perio=le_dom.get_renum_arete_perio();
      // Number of real edges
      int nb_aretes=ok_arete.size();
      int npa=le_dom.numero_premiere_arete();
      for(int i=0; i<nb_aretes; i++)
        {
          if(!ok_arete(i) && tab_pression(npa+i)!=0.)
            {
              Cerr << "Pb: pressure on superfluous edge, P(" << npa+i << ")=" << tab_pression(npa+i) << finl;
              tab_pression(npa+i)=0;
              Process::exit();
            }
          else if ( (renum_arete_perio[i]!=i) && tab_pression(npa+i)!=0.)
            {
              Cerr << "Pb: pressure on superfluous periodic edge, P(" << npa+i << ")=" << tab_pression(npa+i) << finl;
              tab_pression(npa+i)=0;
              Process::exit();
            }
        }
    }
 
  // Apply periodicity on vertices for the pressure:
  if (le_dom.get_alphaS())
    {
      const Domaine& dom=le_dom.domaine();
      int nps=le_dom.numero_premier_sommet();
      int ns=le_dom.domaine().nb_som();
      CIntArrView renum_som_perio = dom.get_renum_som_perio().view_ro();
      DoubleArrView pression = static_cast<DoubleVect&>(tab_pression).view_rw();
      Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                           Kokkos::RangePolicy<>(0, ns), KOKKOS_LAMBDA(
                             const int i)
      {
        int k=renum_som_perio(i);
        if (k!=i) pression(nps+i)=pression(nps+k);
      });
      end_gpu_timer(__KERNEL_NAME__);
    }
  // pression.echange_espace_virtuel();
  // to remove the pressure minimum if no Pref and if only P0, otherwise filter later
  if (le_dom.get_alphaE() && (le_dom.get_alphaS()==0) && (le_dom.get_alphaA()==0) )
    Assembleur_P_VEF::modifier_solution(tab_pression);
  // Verification possible via environment variable:
  char* theValue = getenv("TRUST_VERIFIE_DIRICHLET");
  if(theValue != nullptr) verifier_dirichlet();
 
  return 1;
}
 
void Assembleur_P_VEFPreP1B::verifier_dirichlet()
{
  if (domaine_Vef().get_alphaE()+domaine_Vef().get_alphaS()+domaine_Vef().get_alphaA()!=dimension)
    {
      Cerr << "Assembleur_P_VEFPreP1B::verifier_dirichlet does not yet work with your discretization" << finl;
      Process::exit();
    }
  // Various verifications of boundary conditions
  // by post-processing the result in the Divergence_U field
  IntVect Faces_de_Dirichlet(domaine_Vef().nb_elem_tot());
  IntTab faces(domaine_Vef().nb_elem_tot(),2);
  faces=-1;
  const IntTab& face_sommets=domaine_Vef().face_sommets();
  Faces_de_Dirichlet=0;
  const Conds_lim& les_cl = le_dom_Cl_VEF->les_conditions_limites();
  for(int i=0; i<les_cl.size(); i++)
    {
      const Cond_lim& la_cl = les_cl[i];
      const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
      int num1 = le_bord.num_premiere_face();
      int num2 = num1 + le_bord.nb_faces();
      if ((sub_type(Dirichlet, la_cl.valeur()))
          || (sub_type(Dirichlet_homogene, la_cl.valeur())))
        {
          for (int face=num1; face<num2; face++)
            {
              int elem=domaine_Vef().face_voisins(face,0);
              if (elem!=-1)
                {
                  Faces_de_Dirichlet(elem)++;
                  if (faces(elem,0)==-1) faces(elem,0)=face;
                  else faces(elem,1)=face;
 
                }
              elem=domaine_Vef().face_voisins(face,1);
              if (elem!=-1) Faces_de_Dirichlet(elem)++;
            }
        }
    }
  DoubleVect& post=ref_cast(Navier_Stokes_std,mon_equation.valeur()).div().valeurs();
  post=0;
  for (int elem=0; elem<Faces_de_Dirichlet.size_array(); elem++)
    {
      //post(elem)=Faces_de_Dirichlet(elem);
      if (Faces_de_Dirichlet(elem)>1)
        {
          Cerr << "Element " << elem << " has " << Faces_de_Dirichlet(elem) << " Dirichlet faces separated by a ";
          const IntVect& ok_arete=domaine_Vef().get_ok_arete();
          const IntTab& aretes_som=domaine_Vef().domaine().aretes_som();
          // Loop over edges to verify boundary edges
          for (int k=0; k<6; k++)
            {
              int arete=domaine_Vef().domaine().elem_aretes(elem,k);
              // The 2 vertices of the edge
              int S0=aretes_som(arete,0);
              int S1=aretes_som(arete,1);
              // Verify if they appear twice in the vertices of the faces,
              // in which case this is the edge shared by 2 Dirichlet faces
              int ok=0;
              for (int l=0; l<2; l++)
                {
                  int face=faces(elem,l);
                  for (int s=0; s<3; s++)
                    {
                      if (face_sommets(face,s)==S0) ok++;
                      if (face_sommets(face,s)==S1) ok++;
                    }
                }
              if (ok==4) Cerr << "edge " << (ok_arete(arete)==0?"superfluous":"") << finl;
              // Count the superfluous edges
              post(elem)+=(ok_arete(arete)==0);
            }
          //Process::exit();
        }
    }
}
 
void Assembleur_P_VEFPreP1B::projete_L2(DoubleTab& pression)
{
  if (domaine_Vef().get_alphaE()+domaine_Vef().get_alphaS()+domaine_Vef().get_alphaA()!=3)
    {
      Cerr << "Assembleur_P_VEFPreP1B::projete_L2 only works in P0+P1+Pa mode" << finl;
      Process::exit();
    }
  //Cerr << "Projection L2" << finl;
  const Domaine_VEF& domaine_vef = domaine_Vef();
  const Domaine& domaine = domaine_vef.domaine();
  int nb_elem_tot=domaine.nb_elem_tot();
  int ns=domaine.nb_som_tot();
  int nb_elem=domaine.nb_elem();
  int i,j;
  const IntTab& som_elem=domaine_vef.domaine().les_elems();
  const DoubleVect& volumes=domaine_vef.volumes();
  double volume_tot=0;
  double somme=0;
  int nsr=nb_elem_tot+ns;
  int nse=domaine.nb_som_elem();
  for(i=nb_elem_tot; i<nsr; i++)
    somme+=pression(i);
  somme/=ns;
  for(i=nb_elem_tot; i<nsr; i++)
    pression(i)-=somme;
  static double dalpha =(1./dimension*(dimension+1));
  double pmoy=0;
  for(i=0; i<nb_elem; i++)
    {
      double p=pression(i);
      double v=volumes(i);
      for(j=0; j<nse; j++)
        p+=pression(nb_elem_tot+som_elem(i,j))*dalpha;
      pmoy+=p*v;
      volume_tot+=v;
    }
  pmoy/=volume_tot;
  for(i=0; i<nb_elem; i++)
    pression(i)-=pmoy;
}
 
/*! @brief Optionally modify the matrix to make it definite if it is not. Default values:
 *
 *     Constraints:
 *     Access: input
 *
 * @return 1 if the matrix was modified, 0 otherwise.
 */
int Assembleur_P_VEFPreP1B::modifier_matrice(Matrice& la_matrice)
{
  int matrice_modifiee=0;
  has_P_ref=0;
  const Conds_lim& les_cl = le_dom_Cl_VEF->les_conditions_limites();
  const Domaine_VEF& domaine_VEF = domaine_Vef();
  // Check if there is a reference pressure; if so, the matrix is not modified
  for(int i=0; i<les_cl.size(); i++)
    if (sub_type(Neumann_sortie_libre,les_cl[i].valeur()) ||sub_type(Robin_VEF, les_cl[i].valeur()))
      {
        has_P_ref=1;
        if (domaine_VEF.get_alphaA())
          {
            // Also verify that the pressure is indeed zero on Pa support
            const DoubleTab& val=ref_cast(Neumann_sortie_libre,les_cl[i].valeur()).champ_front().valeurs();
            int nbval=val.dimension(0);
            for (int n=0; n<nbval; n++)
              if (val(n,0)!=0)
                {
                  Cerr << "Non-zero imposed pressure boundary condition is not yet" << finl;
                  Cerr << "supported in VEF with edge support Pa." << finl;
                  // Work to be done in Assembleur_P_VEFPreP1B::modifier_secmem_aretes
                  Process::exit();
                }
          }
      }
 
  Matrice_Bloc_Sym& matrice=ref_cast(Matrice_Bloc_Sym, la_matrice.valeur());
  int P0 = 0;
  int P1 = P0 + domaine_VEF.get_alphaE();
  int Pa = P1 + domaine_VEF.get_alphaS();
 
 
 
  int CL_neumann=has_P_ref;
  ////////////
  // Partie P0
  ////////////
  if (domaine_VEF.get_alphaE())
    {
      // Impose a reference pressure on an element if no Neumann boundary condition
      if (has_P_ref)
        {
          Matrice_Bloc& mat_bloc=ref_cast(Matrice_Bloc,matrice.get_bloc(P0,P0).valeur());
          Matrice_Morse_Sym& A00RR = ref_cast(Matrice_Morse_Sym,mat_bloc.get_bloc(0,0).valeur());
          A00RR.set_est_definie(1);
        }
      matrice_modifiee=Assembleur_P_VEF::modifier_matrice(matrice.get_bloc(P0,P0));
    }
 
  ////////////
  // Partie P1
  ////////////
  if (domaine_VEF.get_alphaS())
    {
      Matrice_Bloc& mat_bloc_p1_p1 = ref_cast(Matrice_Bloc, matrice.get_bloc(P1,P1).valeur());
      Matrice_Morse_Sym& A11RR = ref_cast(Matrice_Morse_Sym,mat_bloc_p1_p1.get_bloc(0,0).valeur());
      // Impose a reference pressure on a vertex if P0 support or if no Neumann boundary condition
      const bool is_first_proc_with_real_elems = Process::me() == Process::mp_min(le_dom_VEF->nb_elem() ? Process::me() : 1e8);
      if (((!(Process::mp_max(CL_neumann))) || ((domaine_VEF.get_alphaE()) && (domaine_VEF.get_cl_pression_sommet_faible()==1) ) ) && is_first_proc_with_real_elems)
        {
          int sommet_referent=0;
          double distance=DMAXFLOAT;
          const DoubleTab& coord=domaine_VEF.domaine().coord_sommets();
          int nb_som=domaine_VEF.domaine().nb_som();
          for(int i=0; i<nb_som; 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(A11RR(i,i),0.))
                {
                  distance=tmp;
                  sommet_referent=i;
                }
            }
          //Cerr << "Modifying vertex row " << sommet_referent << finl;
          A11RR(sommet_referent,sommet_referent)*=2;
          // has_P_ref=1;
          matrice_modifiee=1;
        }
      A11RR.set_est_definie(1);
    }
 
  ////////////
  // Partie Pa
  ////////////
  if (domaine_VEF.get_alphaA())
    {
      Matrice_Bloc& mat_bloc_pa_pa = ref_cast(Matrice_Bloc, matrice.get_bloc(Pa,Pa).valeur());
      Matrice_Morse_Sym& A22RR = ref_cast(Matrice_Morse_Sym,mat_bloc_pa_pa.get_bloc(0,0).valeur());
      // Impose a reference pressure on an edge in P0+Pa mode only
      const bool is_first_proc_with_real_elems = Process::me() == Process::mp_min(le_dom_VEF->nb_elem() ? Process::me() : 1e8);
      if ((domaine_VEF.get_alphaE() && !domaine_VEF.get_alphaS()) && is_first_proc_with_real_elems)
        {
          int arete_referente=0;
          double distance=DMAXFLOAT;
          const DoubleTab& coord=domaine_VEF.xa();
          int nb_aretes=domaine_VEF.domaine().nb_aretes();
          for(int i=0; i<nb_aretes; i++)
            {
              double tmp=0;
              for (int j=0; j<dimension; j++)
                tmp+=coord(i,j)*coord(i,j);
              if (tmp<distance && A22RR(i,i)!=0)
                {
                  distance=tmp;
                  arete_referente=i;
                }
            }
          //Cerr << "Modifying edge row " << arete_referente << finl;
          A22RR(arete_referente,arete_referente)*=2;
          //has_P_ref=1;
          matrice_modifiee=1;
        }
      A22RR.set_est_definie(1);
    }
 
  {
    Matrice_Bloc& bloc_P0_P0 = ref_cast( Matrice_Bloc,matrice.get_bloc(P0,P0).valeur( ) );
    Matrice_Morse_Sym& A00RR = ref_cast( Matrice_Morse_Sym, bloc_P0_P0.get_bloc(0,0).valeur() );
    matrice.set_est_definie( A00RR.get_est_definie( ) );
  }
 
  return matrice_modifiee;
}
 
inline void range(double& prod, int& i, int& n, int& j, int& m, Matrice_Morse& ARR, Matrice_Morse& ARV, Matrice_Morse& AVR, Matrice_Morse& AVV)
{
  if (i<n)
    if (j<m)
      ARR(i,j)+=prod;
    else
      ARV(i,j-m)+=prod;
  else if (j<m)
    AVR(i-n,j)+=prod;
  else
    AVV(i-n,j-m)+=prod;
}
 
void operation11(Matrice_Bloc& A00, Matrice_Bloc& A01, Matrice_Bloc& A11, double beta, const Domaine& domaine)
{
  //Cerr << "Operation11" << finl;
  Matrice_Morse_Sym& A11RR=ref_cast(Matrice_Morse_Sym, A11.get_bloc(0,0).valeur());
  Matrice_Morse& A11RV=ref_cast(Matrice_Morse, A11.get_bloc(0,1).valeur());
  Matrice_Morse& A11VR=ref_cast(Matrice_Morse, A11.get_bloc(1,0).valeur());
  Matrice_Morse& A11VV=ref_cast(Matrice_Morse, A11.get_bloc(1,1).valeur());
  const IntTab& les_elems=domaine.les_elems();
  const Domaine& dom=domaine;
  int nb_som=A11RR.nb_lignes();
  int nb_som_elem=les_elems.dimension(1);
  // Traverse the elements of matrix A00
  //Cerr << "[" << Process::me() << "] Contribution de A00 dans A11~" << finl;
  int ligne=0;
  for (int i_bloc=0; i_bloc<A00.nb_bloc_lignes(); i_bloc++)
    {
      int colonne=0;
      int nb_lignes=0;
      for (int j_bloc=0; j_bloc<A00.nb_bloc_colonnes(); j_bloc++)
        {
          Matrice_Morse& A00ij=ref_cast(Matrice_Morse, A00.get_bloc(i_bloc,j_bloc).valeur());
          const auto* tab1=A00ij.get_tab1().addr();
          const int* tab2=A00ij.get_tab2().addr();
          const double* coeff=A00ij.get_coeff().addr();
          nb_lignes=A00ij.nb_lignes();
          int nb_colonnes=A00ij.nb_colonnes();
          for (int i=0; i<nb_lignes; i++)
            {
              int k1=ligne+i; // Element k1
              for (auto n=tab1[i]-1; n<tab1[i+1]-1; n++)
                {
                  int k2=colonne+tab2[n]-1; // Element k2
                  if (k2>=k1)
                    {
                      double prod = beta * beta * coeff[n];        // Compute beta*beta*Ak1k2
                      for (int som1=0; som1<nb_som_elem; som1++)
                        {
                          int s1 = dom.get_renum_som_perio(les_elems(k1,som1));
                          for (int som2=0; som2<nb_som_elem; som2++)
                            {
                              int s2 = dom.get_renum_som_perio(les_elems(k2,som2));
                              if (s2>=s1) range(prod,s1,nb_som,s2,nb_som,A11RR,A11RV,A11VR,A11VV);
                            }
                        }
                      if (k1!=k2)
                        {
                          for (int som1=0; som1<nb_som_elem; som1++)
                            {
                              int s1 = dom.get_renum_som_perio(les_elems(k2,som1));
                              for (int som2=0; som2<nb_som_elem; som2++)
                                {
                                  int s2 = dom.get_renum_som_perio(les_elems(k1,som2));
                                  if (s2>=s1) range(prod,s1,nb_som,s2,nb_som,A11RR,A11RV,A11VR,A11VV);
                                }
                            }
                        }
                    }
                }
            }
          colonne+=nb_colonnes;
        }
      ligne+=nb_lignes;
    }
  // Traverse the elements of matrix A01
  //Cerr << "[" << Process::me() << "] Contribution de A01 dans A11~" << finl;
  ligne=0;
  for (int i_bloc=0; i_bloc<A01.nb_bloc_lignes(); i_bloc++)
    {
      int colonne=0;
      int nb_lignes=0;
      for (int j_bloc=0; j_bloc<A01.nb_bloc_colonnes(); j_bloc++)
        {
          Matrice_Morse& A01ij=ref_cast(Matrice_Morse, A01.get_bloc(i_bloc,j_bloc).valeur());
          const auto* tab1=A01ij.get_tab1().addr();
          const int* tab2=A01ij.get_tab2().addr();
          const double* coeff=A01ij.get_coeff().addr();
          nb_lignes=A01ij.nb_lignes();
          int nb_colonnes=A01ij.nb_colonnes();
          for (int i=0; i<nb_lignes; i++)
            {
              int k=ligne+i; // Element k
              for (auto n=tab1[i]-1; n<tab1[i+1]-1; n++)
                {
                  int s1 = dom.get_renum_som_perio(colonne+tab2[n]-1); // Vertex s1
                  double prod = -beta * coeff[n];        // Compute -beta*Aks
                  for (int som=0; som<nb_som_elem; som++)
                    {
                      int s2 = dom.get_renum_som_perio(les_elems(k,som)); // Vertex s2
                      if (s2>=s1) range(prod,s1,nb_som,s2,nb_som,A11RR,A11RV,A11VR,A11VV);
                      if (s1>=s2) range(prod,s2,nb_som,s1,nb_som,A11RR,A11RV,A11VR,A11VV);
                    }
                }
            }
          colonne+=nb_colonnes;
        }
      ligne+=nb_lignes;
    }
}
 
void operation01(Matrice_Bloc& A00, Matrice_Bloc& A01, double alpha, double beta, const Domaine& domaine)
{
  //Cerr << "[" << Process::me() << "] Operation01" << finl;
  Matrice_Morse& A01RR=ref_cast(Matrice_Morse, A01.get_bloc(0,0).valeur());
  Matrice_Morse& A01RV=ref_cast(Matrice_Morse, A01.get_bloc(0,1).valeur());
  Matrice_Morse& A01VR=ref_cast(Matrice_Morse, A01.get_bloc(1,0).valeur());
  Matrice_Morse& A01VV=ref_cast(Matrice_Morse, A01.get_bloc(1,1).valeur());
  const IntTab& les_elems=domaine.les_elems();
  const Domaine& dom=domaine;
  int nb_elem=A01RR.nb_lignes();
  int nb_som=A01RR.nb_colonnes();
  int nb_som_elem=les_elems.dimension(1);
  // Iterate over the coefficients of A00
  int ligne=0;
  for (int i_bloc=0; i_bloc<A00.nb_bloc_lignes(); i_bloc++)
    {
      int colonne=0;
      int nb_lignes=0;
      for (int j_bloc=0; j_bloc<A00.nb_bloc_colonnes(); j_bloc++)
        {
          Matrice_Morse& A00ij=ref_cast(Matrice_Morse, A00.get_bloc(i_bloc,j_bloc).valeur());
          const auto* tab1=A00ij.get_tab1().addr();
          const int* tab2=A00ij.get_tab2().addr();
          const double* coeff=A00ij.get_coeff().addr();
          nb_lignes=A00ij.nb_lignes();
          int nb_colonnes=A00ij.nb_colonnes();
          int s;
          for (int i=0; i<nb_lignes; i++)
            {
              int k1=ligne+i; // Element k1
              for (auto n=tab1[i]-1; n<tab1[i+1]-1; n++)
                {
                  int k2=colonne+tab2[n]-1; // Element k2
                  if (k2>=k1)
                    {
                      double prod = -alpha * beta * coeff[n];        // Compute -alpha*beta*Ak1k2
                      for (int som=0; som<nb_som_elem; som++)
                        {
                          s = dom.get_renum_som_perio(les_elems(k2,som));
                          range(prod,k1,nb_elem,s,nb_som,A01RR,A01RV,A01VR,A01VV);
                          if (k1!=k2)
                            {
                              s = dom.get_renum_som_perio(les_elems(k1,som));
                              range(prod,k2,nb_elem,s,nb_som,A01RR,A01RV,A01VR,A01VV);
                            }
                        }
                    }
                }
            }
          colonne+=nb_colonnes;
        }
      ligne+=nb_lignes;
    }
}
 
void Assembleur_P_VEFPreP1B::changer_base_matrice(Matrice& la_matrice)
{
  assert(domaine_Vef().get_alphaE() && domaine_Vef().get_alphaS() && !domaine_Vef().get_alphaA()); // P0+P1 uniquement
  Cerr << "Basis change for the matrix: P0+P1->P1Bulle" << finl;
  Matrice_Bloc_Sym& matrice=ref_cast(Matrice_Bloc_Sym, la_matrice.valeur());
  Matrice_Bloc& A00=ref_cast(Matrice_Bloc, matrice.get_bloc(0,0).valeur());
  Matrice_Bloc& A01=ref_cast(Matrice_Bloc, matrice.get_bloc(0,1).valeur());
  Matrice_Bloc& A11=ref_cast(Matrice_Bloc, matrice.get_bloc(1,1).valeur());
 
  // Modification of block A11
  // As1s2~=As1s2-beta*[somme(Ak1s2)(s1 belongs to k1)+somme(Ak1s1)(s2 belonging to k1)]+beta*beta*somme(Ak1k2)(s1 belonging to k1 and s2 belonging to k2)
  operation11(A00,A01,A11,beta_,domaine_Vef().domaine());
 
  // Modification of block A01
  // Ak1s~=alpha*Ak1s - alpha*beta*somme(Ak1k2)(s belonging to k2)
  A01*=alpha_;
  operation01(A00,A01,alpha_,beta_,domaine_Vef().domaine());
 
  // Modification of block A00
  // Ak1k2~ = alpha * alpha * Ak1k2
  A00*=alpha_*alpha_;
}
 
 
void Assembleur_P_VEFPreP1B::changer_base_second_membre(DoubleVect& y)
{
  // ys~ = ys - beta * somme(yk)(s belonging to k)
  // yk~ = alpha * yk
  changer_base<vecteur::second_membre>(y);
  y.echange_espace_virtuel();
}
 
void Assembleur_P_VEFPreP1B::changer_base_pression_inverse(DoubleVect& x)
{
  // xk = alpha * xk~ - beta * somme(xs~)(s belonging to k)
  // xs = xs~
  changer_base<vecteur::pression_inverse>(x);
  assert(check_espace_virtuel_vect(x));
}
 
void Assembleur_P_VEFPreP1B::changer_base_pression(DoubleVect& x)
{
  // xk~ = xk / alpha + beta / alpha * somme(xs)(s belonging to k)
  // xs~ = xs
  changer_base<vecteur::pression>(x);
  assert(check_espace_virtuel_vect(x));
}
 
template<vecteur _v_>
void Assembleur_P_VEFPreP1B::changer_base(DoubleVect& tab_v)
{
  assert(domaine_Vef().get_alphaE() && domaine_Vef().get_alphaS() && !domaine_Vef().get_alphaA()); // P0+P1 uniquement
  // xk~ = xk / alpha + beta / alpha * somme(xs)(s belonging to k)
  // xs~ = xs
  double alpha = alpha_;
  double beta = beta_;
  int nb_elem_tot = domaine_Vef().nb_elem_tot();
  int nb_som_elem = domaine_Vef().domaine().les_elems().dimension(1);
  CIntTabView les_elems = domaine_Vef().domaine().les_elems().view_ro();
  CIntArrView renum_som_perio = domaine_Vef().domaine().get_renum_som_perio().view_ro();
  DoubleArrView v = tab_v.view_rw();
  Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nb_elem_tot, KOKKOS_LAMBDA(const int k)
  {
    if (_v_==vecteur::pression) v(k) /= alpha;
    double somme=0;
    for (int som=0; som<nb_som_elem; som++)
      {
        int s = nb_elem_tot + renum_som_perio(les_elems(k, som));
        if (_v_==vecteur::pression)         somme += v(s);
        if (_v_==vecteur::pression_inverse) somme += v(s);
        if (_v_==vecteur::second_membre) Kokkos::atomic_sub(&v(s), beta * v(k));
      }
    if (_v_==vecteur::pression)         v(k) += beta / alpha * somme;
    if (_v_==vecteur::pression_inverse) v(k) = alpha * v(k) - beta * somme;
    if (_v_==vecteur::second_membre)    v(k) *= alpha;
  });
  end_gpu_timer(__KERNEL_NAME__);
}