Champ_P1iP1B_implementation.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 <Champ_P1iP1B_implementation.h>
#include <Champ_implementation_P1.h>
#include <Domaine.h>
#include <Domaine_VEF.h>
#include <Domaine_Cl_VEF.h>
#include <Debog.h>
 
#include <Matrice_Morse_Sym.h>
#include <TRUSTLists.h>
#include <SolveurSys.h>
#include <Solv_GCP.h>
#include <SSOR.h>
 
#include <Neumann_sortie_libre.h>
#include <Champ_P0_VEF.h>
#include <TRUSTTab_parts.h>
#include <Frontiere_dis_base.h>
 
#include <kokkos++.h>
#include <TRUSTArray_kokkos.tpp>
 
DoubleVect& Champ_P1iP1B_implementation::valeur_a_elem(const DoubleVect& position, DoubleVect& val, int le_poly) const
{
  int dimension=Objet_U::dimension;
  DoubleTab vals(1,1);
  IntVect les_polys(1);
  DoubleTab positions(1,dimension);
 
  les_polys(0) = le_poly;
  for (int i=0; i<dimension; i++)
    positions(0,i) = position(i);
 
  valeur_aux_elems(positions, les_polys, vals);
 
  val(0) = vals(0,0);
  return val;
}
 
double Champ_P1iP1B_implementation::valeur_a_elem_compo(const DoubleVect& position, int le_poly, int ncomp) const
{
  double val=0;
  int dimension=Objet_U::dimension;
 
  DoubleVect vals(1);
  IntVect les_polys(1);
  DoubleTab positions(1,dimension);
 
  les_polys(0) = le_poly;
  for (int i=0; i<dimension; i++)
    positions(0,i) = position(i);
 
  valeur_aux_elems_compo(positions, les_polys, vals, ncomp);
 
  val = vals(0);
 
  return val;
}
 
// Recupere un tableau positions contenant des coordonnees de points
// contenus dans les elements ranges dans le tableau les_polys
// Renvoie le tableau val contenant les valeurs du champ aux points
DoubleTab& Champ_P1iP1B_implementation::valeur_aux_elems(const DoubleTab& positions, const IntVect& les_polys, DoubleTab& val) const
{
  // zs bien initialise en 3D
 
  const Domaine_VEF& zvef = domaine_vef();
  const Domaine& domaine_geom = zvef.domaine();
  const DoubleTab& coord = domaine_geom.coord_sommets();
  const IntTab& sommet_poly = domaine_geom.les_elems();
  int prs=zvef.numero_premier_sommet();
  int nb_compo_=le_champ().nb_comp();
  int dimension=Objet_U::dimension;
  const int les_polys_size=les_polys.size();
 
  // TODO : FIXME
  // For FT the resize should be done in its good position and not here ...
  if (val.nb_dim() == 1 ) val.resize(les_polys_size,1);
 
  assert((val.dimension(0) == les_polys.size())||(val.dimension_tot(0) == les_polys.size()));
  assert(val.line_size() == nb_compo_);
 
  // Filtrage du champ dans le tableau champ_filtre_
  champ_filtre_=filtrage(zvef,le_champ());
 
  if (nb_compo_ == 1)
    {
      // HIP Kokkos prevent from the function calls inside the Kokkos region:
      // Error: error: reference to __host__ function 'get_alphaE' in __host__ __device__ function
      bool alphaE = zvef.get_alphaE();
      bool alphaS = zvef.get_alphaS();
      CDoubleArrView champ_filtre_v = static_cast<DoubleVect&>(champ_filtre_).view_ro();
      CIntTabView sommet_poly_v = sommet_poly.view_ro();
      CDoubleTabView coord_v = coord.view_ro();
      CDoubleTabView positions_v = positions.view_ro();
      CIntArrView les_polys_v = les_polys.view_ro();
      DoubleTabView val_v = val.view_rw();
      Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),les_polys_size, KOKKOS_LAMBDA(
                             const int rang_poly)
      {
        // On initialise
        val_v(rang_poly,0) = 0.;
 
        // On prend le_poly
        int le_poly=les_polys_v(rang_poly);
        if (le_poly != -1)
          {
            // Contribution P0
            if (alphaE)
              val_v(rang_poly, 0) += champ_filtre_v(le_poly);
            // Contribution P1
            if (alphaS)
              {
                double xs = positions_v(rang_poly,0);
                double ys = positions_v(rang_poly,1);
                double zs = (dimension == 3) ? positions_v(rang_poly,2) : 0;
 
                // Calcul par les fonctions de forme P1 du champ filtre
                // au point (xs,yz,zs)
                for (int i=0; i<dimension+1; i++)
                  {
                    int som = prs+sommet_poly_v(le_poly,i);
                    double champ_filtre_som = champ_filtre_v(som);
                    double li=0.;
 
                    if (dimension == 2)
                      li=coord_barycentrique_P1_triangle(sommet_poly_v,
                                                         coord_v, xs, ys,
                                                         le_poly, i);
                    else if (dimension == 3)
                      li=coord_barycentrique_P1_tetraedre(sommet_poly_v,
                                                          coord_v, xs, ys, zs,
                                                          le_poly, i);
                    val_v(rang_poly,0) += champ_filtre_som * li;
                  }
              }
          }
      });
      end_gpu_timer(__KERNEL_NAME__);
    }
  else // nb_compo_ > 1
    {
      Cerr << "Vous en etes deja la vous ?" << finl;
      Process::exit();
    }
  return val;
}
 
// Idem que valeur_aux_elems mais pour une composante d'un champ vecteur
DoubleVect& Champ_P1iP1B_implementation::valeur_aux_elems_compo(const DoubleTab& positions, const IntVect& les_polys, DoubleVect& val, int ncomp) const
{
  Cerr << "Champ_P1iP1B_implementation::valeur_aux_elems_compo non code." << finl;
  // Factoriser le code avec la methode valeur_aux_elems()
  Process::exit();
  return val;
}
 
 
// Recupere un domaine
// Renvoie un tableau contenant les valeurs du champ aux sommets du domaine
DoubleTab& Champ_P1iP1B_implementation::valeur_aux_sommets(const Domaine& dom, DoubleTab& tab_val) const
{
  const Domaine_VEF& zvef = domaine_vef();
  int nb_compo_=le_champ().nb_comp();
  assert(nb_compo_ == tab_val.line_size());
 
  // Filtrage du champ pour le postraitement (contenu dans le tableau champ_filtre_)
  champ_filtre_=filtrage(zvef,le_champ());
 
  if (zvef.get_alphaE()) // Support P0
    {
      Champ_P0_VEF tmp;
      tmp.associer_domaine_dis_base(zvef);
      tmp.valeurs() = champ_filtre_;
      tmp.valeur_aux_sommets(zvef.domaine(), tab_val);
    }
 
  if (zvef.get_alphaS()) // Support P1
    {
      int nbs=zvef.nb_som();
      int prs=zvef.numero_premier_sommet();
      if (nb_compo_ == 1)
        {
          CDoubleArrView champ_filtre = static_cast<const DoubleVect&>(champ_filtre_).view_ro();
          DoubleTabView val = tab_val.view_rw();
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), nbs, KOKKOS_LAMBDA(const int num_som)
          {
            int som = prs+num_som;
            val(num_som,0) += champ_filtre(som);
          });
          end_gpu_timer(__KERNEL_NAME__);
        }
      else // nb_compo_ > 1
        {
          Cerr << "Vous en etes deja la vous ?" << finl;
          Process::exit();
        }
    }
  return tab_val;
}
 
 
//
DoubleVect& Champ_P1iP1B_implementation::
valeur_aux_sommets_compo(const Domaine& dom,
                         DoubleVect& val,
                         int ncomp) const
{
  Cerr << "Champ_P1iP1B_implementation::valeur_aux_sommets_compo pas code." << finl;
  // Factoriser avec la methode valeur_aux_sommets()
  Process::exit();
  return val;
}
 
DoubleTab& Champ_P1iP1B_implementation::remplir_coord_noeuds(DoubleTab& coord) const
{
  const Domaine_VEF& zvef=domaine_vef();
  const Domaine& le_dom=zvef.domaine();
  const Domaine& dom=le_dom;
  const DoubleTab& coord_sommets=dom.coord_sommets();
  int dimension=Objet_U::dimension;
 
  const DoubleTab& xg = zvef.xp();
  int nbe=zvef.get_alphaE() ? zvef.nb_elem_tot() : 0;
  int nbs=zvef.get_alphaS() ? zvef.nb_som_tot() : 0;
  coord.resize(nbe+nbs,dimension);
  int k = 0;
  if (zvef.get_alphaE())
    {
      for(int i=0; i<nbe; i++)
        {
          for(int j=0; j<dimension; j++)
            coord(k,j)=xg(i,j);
          k++;
        }
    }
  if (zvef.get_alphaS())
    {
      for(int i=0; i<nbs; i++)
        {
          for(int j=0; j<dimension; j++)
            coord(k,j)=coord_sommets(i,j);
          k++;
        }
    }
  return coord;
}
 
int Champ_P1iP1B_implementation::
remplir_coord_noeuds_et_polys(DoubleTab& positions,
                              IntVect& polys) const
{
  Cerr << "Pas code " << finl;
  assert(0);
  Process::exit();
  return 1;
}
 
void assembler(const Domaine_VEF& domaine_VEF, Matrice& matrice)
{
  matrice.typer("Matrice_Morse_Sym");
  Matrice_Morse_Sym& MatPoisson=ref_cast(Matrice_Morse_Sym, matrice.valeur());
  int nb_som_tot = domaine_VEF.domaine().nb_som_tot();
  int nb_arete_tot = domaine_VEF.domaine().nb_aretes_tot();
  int nnz=nb_som_tot+nb_arete_tot;
  const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
  const ArrOfInt& renum_arete_perio=domaine_VEF.get_renum_arete_perio();
  const Domaine& dom=domaine_VEF.domaine();
 
  IntLists voisins(nb_som_tot);
  DoubleLists coeffs(nb_som_tot);
  DoubleVect diag(nb_som_tot);
  // On parcourt toutes les aretes non periodiques:
  for(int arete=0; arete<nb_arete_tot; arete++)
    {
      if (renum_arete_perio[arete]==arete)
        {
          int som1=dom.get_renum_som_perio(aretes_som(arete, 0));
          int som2=dom.get_renum_som_perio(aretes_som(arete, 1));
          if(som1>som2)
            {
              int tmp=som1;
              som1=som2;
              som2=tmp;
            }
          voisins[som1].add(som2);
          coeffs[som1].add(1);
          diag(som1)++;
          diag(som2)++;
        }
    }
  for (int i=0; i<nb_som_tot; i++)
    {
      if (diag(i)==0)
        {
          assert(i!=dom.get_renum_som_perio(i));
          diag(i)=1; // Sommets periodiques
        }
    }
  MatPoisson.dimensionner(nb_som_tot, nnz) ;
  MatPoisson.remplir(voisins, coeffs, diag) ;
  MatPoisson.compacte() ;
  MatPoisson.set_est_definie(1);
  //Cerr << "Matrice P1";MatPoisson.imprimer(Cerr);
}
 
static double coeff = 4. / 5.;
static double coeff_inv = 1 / coeff;
double second_membre(const Domaine_VEF& domaine_VEF, ArrOfDouble& Pa, DoubleVect& secmem)
{
  int nb_arete_tot = domaine_VEF.domaine().nb_aretes_tot();
  const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
  const ArrOfInt& renum_arete_perio=domaine_VEF.get_renum_arete_perio();
  const Domaine& dom=domaine_VEF.domaine();
  secmem=0;
  // On parcourt toutes les aretes non periodiques
  for(int arete=0; arete<nb_arete_tot; arete++)
    {
      if (renum_arete_perio[arete]==arete)
        {
          double x=Pa[arete];
          int som1=dom.get_renum_som_perio(aretes_som(arete, 0));
          int som2=dom.get_renum_som_perio(aretes_som(arete, 1));
          secmem(som1)+=x;
          secmem(som2)+=x;
        }
    }
  secmem*=coeff;
  secmem.echange_espace_virtuel();
  double norme=mp_norme_vect(secmem);
  //Cerr << "||somme Pa||= " << norme << finl;
  return norme;
}
 
void corriger(const Domaine_VEF& domaine_VEF, DoubleTab& champ_filtre_, Matrice& matrice, const int Condition_Neumann_imposee_)
{
  int nb_elem=domaine_VEF.nb_elem();
  int nb_som=domaine_VEF.nb_som();
 
  // Tableaux d'acces aux valeurs Pk et Ps
  if (!domaine_VEF.get_alphaE() || !domaine_VEF.get_alphaS())
    {
      Cerr << "Erreur dans Champ_P1iP1B_impl.cpp: corriger(...), le champ doit avoir une partie sommets et elements" << finl;
      Process::exit();
    }
  DoubleTab_parts parties_P(champ_filtre_);
  DoubleVect& Pk = parties_P[0];  // partie elements
  DoubleVect& Ps = parties_P[1];  // partie sommets
 
  // Filtrage si support arete
  if (domaine_VEF.get_alphaA())
    {
      if (domaine_VEF.get_renum_arete_perio().size_array()==0)
        {
          Cerr<<"We try to Champ_P1iP1B_impl::corriger but get_renum_arete_perio"<<finl;
        }
    }
  if (domaine_VEF.get_alphaA()&& domaine_VEF.get_renum_arete_perio().size_array())
    {
      champ_filtre_.ensureDataOnHost(); // We need copy H2D as the Pa part of DoubleTab_parts is computed on host
      ToDo_Kokkos("critical for P0P1Pa");
      DoubleVect& Pa = parties_P[2];  // partie aretes
 
      // Si premier passage on assemble la matrice
      if (!matrice) assembler(domaine_VEF, matrice);
 
      // Construction du second membre
      DoubleVect secmem;
      domaine_VEF.domaine().creer_tableau_sommets(secmem);
      second_membre(domaine_VEF, Pa, secmem);
 
      // Calcul de la correction
      DoubleVect solution;
      solution.copy(secmem, RESIZE_OPTIONS::NOCOPY_NOINIT); // Copie de la structure
      solution = 0.;
 
      SolveurSys solveur;
      solveur.typer("Solv_GCP");
      OWN_PTR(Precond_base) p;
      p.typer("SSOR");
      ref_cast(Solv_GCP,solveur.valeur()).set_precond(p);
      solveur.nommer("Pa_filter_solver");
      solveur.resoudre_systeme(matrice, secmem, solution);
 
      // Application de la periodicite sur la solution:
      const Domaine& dom=domaine_VEF.domaine();
      for(int i=0; i<nb_som; i++)
        {
          int som=dom.get_renum_som_perio(i);
          if (som!=i)
            solution(i)=solution(som);
        }
 
      // Correction de Ps:
      for(int som=0; som<nb_som; som++)
        Ps(som) += solution(som);
 
      // Correction de Pa:
      int nb_arete = domaine_VEF.domaine().nb_aretes();
      const IntTab& aretes_som=domaine_VEF.domaine().aretes_som();
      const ArrOfInt& ok_arete=domaine_VEF.get_ok_arete();
      for(int arete=0; arete<nb_arete; arete++)
        {
          if(!ok_arete[arete] && Pa(arete)!=0)
            {
              Cerr << "Pa(arete_superflue)!=0 dans Champ_P1iP1B_implementation::corriger" << finl;
              Cerr << "Contacter le support TRUST." << finl;
              Cerr << "S'il s'agit d'une reprise au format xyz d'un calcul" << finl;
              Cerr << "alors il faut ecraser le fichier .ok_arete de votre calcul de reprise" << finl;
              Cerr << "par le fichier .ok_arete de votre calcul precedant afin d'avoir les" << finl;
              Cerr << "memes aretes superflues." << finl;
              Process::exit();
            }
          int som1=aretes_som(arete, 0);
          int som2=aretes_som(arete, 1);
          double sigma=solution(som1)+solution(som2);
          Pa(arete)-=sigma*coeff_inv;
        }
 
      // Correction de Pk:
      const IntTab& les_elems = domaine_VEF.domaine().les_elems();
      int nb_som_par_elem=les_elems.dimension_tot(1);
      for(int elem=0; elem<nb_elem; elem++)
        {
          double sigma=0;
          for (int som=0; som<nb_som_par_elem; som++)
            sigma+=solution(les_elems(elem,som));
          Pk(elem)-=sigma/4;
        }
    }
 
  // Filtrage L2:
  double moyenne_K = mp_moyenne_vect(Pk);                // Calcul de la moyenne du champ aux elements
  Pk -= moyenne_K;                                 // Correction des elements du champ_filtre_
  Ps += moyenne_K;                                 // Correction des sommets du champ_filtre_
 
  // On retranche a Ps sa moyenne si pas de Cl de Neumann pour avoir une moyenne nulle
  // Il faudrait peut etre mieux faire cela dans Assembleur_P_VEFPreP1B::modifier_solution
  // si cela est possible...
  assert(Condition_Neumann_imposee_!=-1);
 
  if (!Condition_Neumann_imposee_)
    {
      // L'espace virtuel de Ps n'a pas besoin d'etre a jour pour calculer correctement la moyenne
      Ps -= mp_moyenne_vect(Ps); // Correction des sommets du champ_filtre_
    }
 
  champ_filtre_.echange_espace_virtuel();         // Mise a jour des espaces virtuels
}
 
DoubleTab& Champ_P1iP1B_implementation::filtrage(const Domaine_VEF& zvef, const Champ_base& un_champ, bool implicitCoupling) const
{
  // Filtrage si supports element et sommet presents au moins
  if (zvef.get_alphaE() && zvef.get_alphaS())
    {
      // No filtering if it has already been done on this field during the time step
      // and if no sub-iterations are performed during the time step
      // The condition implicitCoupling=true enables filtering and therefore updates
      // the variable "champ_filtre_" during sub-iterations.
      // Sub-iterations can be performed within a time step in the case of an implicit
      // coupling with another code, for example with a structural solver
      // in fluid-structure interaction applications
      if (un_champ.valeurs().data()==adresse_champ_filtre_ // PL: data() and not addr(). Else copyFromDevice on GPU !
          && un_champ.temps()==temps_filtrage_
          && !implicitCoupling)
        return champ_filtre_;
 
      // On copie le champ a filtrer dans le tableau qui contiendra le champ filtre
      champ_filtre_=un_champ.valeurs();
      temps_filtrage_=un_champ.temps();
      adresse_champ_filtre_=un_champ.valeurs().data();
      //Cout << "Filtrage du champ " << un_champ.le_nom() << " au temps " << un_champ.temps() << finl;
 
      // Correction du champ
      corriger(zvef, champ_filtre_, matrice_filtrage_, Condition_Neumann_imposee_);
      Debog::verifier("champ apres filtre=", champ_filtre_);
      return champ_filtre_;
    }
  else
    {
      champ_filtre_=un_champ.valeurs();
      return champ_filtre_;
    }
}
 
// Fixe Condition_Neumann_imposee_ a 1 si il existe une CL de Neumann 0 sinon
// Si pas de condition de Neumann, on imposera une moyenne nulle dans le filtrage du champ
// pour avoir le meme champ en sequentiel et en parallele
void Champ_P1iP1B_implementation::completer(const Domaine_Cl_dis_base& zcl)
{
  const Conds_lim& les_cl = zcl.les_conditions_limites();
  Condition_Neumann_imposee_ = 0;
  for(int i=0; i<les_cl.size(); i++)
    if (sub_type(Neumann_sortie_libre,les_cl[i].valeur()))
      Condition_Neumann_imposee_ = 1;
}
 
DoubleTab& Champ_P1iP1B_implementation::trace(const Frontiere_dis_base& fr, const DoubleTab& y, DoubleTab& x,int distant) const
{
  assert(distant==0);
  const Front_VF& fr_vf = ref_cast(Front_VF, fr);
  const Domaine_VEF& domaine_dis = domaine_vef();
  const DoubleTab& cg_faces = domaine_dis.xv();
  int nb_faces_fr = fr_vf.nb_faces();
  int face, elem;
 
  assert(x.dimension(0)==fr_vf.nb_faces());
 
  IntVect elem_voisins(nb_faces_fr);
  elem_voisins = -1;
  DoubleTab cg_faces_fr(nb_faces_fr,Objet_U::dimension);
  DoubleTab val_interp(nb_faces_fr, 1);
 
  //On recupere les elements voisins des faces de la frontiere
  for (int i=0; i<nb_faces_fr; i++)
    {
      face = fr_vf.num_premiere_face()+i;
      elem = domaine_dis.face_voisins(face,0);
      if (elem==-1)
        elem = domaine_dis.face_voisins(face,1);
      elem_voisins(i) = elem;
    }
 
  //Remplissage du tableau contenant les coordonnees
  //du centre de gravite des faces de la frontiere
  for (int num_face=0; num_face<nb_faces_fr; num_face++)
    for (int dim=0; dim<Objet_U::dimension; dim++)
      {
        face = fr_vf.num_premiere_face()+num_face;
        cg_faces_fr(num_face,dim) = cg_faces(face,dim);
      }
 
  //Interpolation des valeurs du champ au centre de gravite des faces frontieres
  valeur_aux_elems(cg_faces_fr,elem_voisins,val_interp);
 
  for (int i=0; i<nb_faces_fr; i++)
    {
      x(i,0) = val_interp(i);
    }
  // Useless ? x.echange_espace_virtuel();
  return x;
}