Op_EF_base.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 <Op_EF_base.h>
#include <Domaine_EF.h>
#include <Domaine_Cl_EF.h>
#include <Dirichlet.h>
#include <Dirichlet_homogene.h>
#include <Symetrie.h>
#include <Matrice_Morse.h>
#include <Equation_base.h>
 
#include <Sortie.h>
#include <Operateur.h>
#include <Probleme_base.h>
 
#include <Champ_Uniforme.h>
#include <Schema_Temps_base.h>
#include <Milieu_base.h>
#include <Operateur_base.h>
#include <Operateur_Diff_base.h>
#include <EcrFicPartage.h>
#include <Array_tools.h>
#include <Echange_interne_global_impose.h>
#include <Echange_interne_global_parfait.h>
#include <Champ_front_calc_interne.h>
 
#include <vector>
 
/*! @brief Sizes the matrix that will receive the coefficients from the convection and diffusion operators for the face-based case.
 *
 *  This matrix has a Morse matrix structure.
 *  We start by computing the sizes of the arrays tab1 and tab2.
 *
 */
 
void Op_EF_base::dimensionner(const Domaine_EF& le_dom,
                              const Domaine_Cl_EF& le_dom_cl,
                              Matrice_Morse& la_matrice) const
{
 
  int matrice_stocke=1;
  if (matrice_stocke)
    {
      //if (matrice_sto_.nb_colonnes()) PL: not wrong in //, nb_colonnes() may be zero locally
      if (Process::mp_sum(matrice_sto_.nb_colonnes()))
        {
          la_matrice=matrice_sto_;
          return;
        }
 
    }
  // Sizing the matrix that will receive coefficients from
  // convection and diffusion in the face case.
  // This matrix has a Morse matrix structure.
  // We start by computing the sizes of arrays tab1 and tab2.
  // For this we need to find the neighboring vertices of the considered vertex.
 
  int nb_som=le_dom.nb_som();
  int nfin = le_dom.nb_som_tot();
  // essai
  nb_som=nfin;
 
  const DoubleTab& champ_inconnue = le_dom_cl.equation().inconnue().valeurs();
  int nb_comp = champ_inconnue.line_size();
 
  // for now we don't handle components
 
  const IntTab& elems=le_dom.domaine().les_elems();
  //int nb_elem=le_dom.domaine().nb_elem();
  int nb_elem_tot=le_dom.domaine().nb_elem_tot();
  int nb_som_elem=le_dom.domaine().nb_som_elem();
  int nb_som_face=le_dom.nb_som_face();
 
  // Computation of the number of non void coeffs:
  int extra_nb_coeff = 0;
  // Specific BC that will require extra terms in the matrix:
  for(int i=0; i<le_dom_cl.nb_cond_lim(); i++)
    {
      const Cond_lim_base& ccl = le_dom_cl.les_conditions_limites(i).valeur();
      if(sub_type(Echange_interne_global_impose, ccl))
        {
          if (nb_som_face != 1)
            {
              Cerr << "Boundary condition 'Echange_interne_global_impose' or 'Echange_interne_parfait' is not validated for EF and mesh_dim > 1 (i.e. sth other than a wire mesh)!!" << finl;
              Process::exit(-1);
            }
          const Echange_interne_global_impose& cl = ref_cast(Echange_interne_global_impose, ccl);
          const Champ_front_calc_interne& ch = ref_cast(Champ_front_calc_interne, cl.T_ext());
          extra_nb_coeff += ch.face_map().size_array();
          if(sub_type(Echange_interne_global_parfait, ccl))
            extra_nb_coeff += ch.face_map().size_array() / 2;
        }
    }
  int nb_coeff=(int)(nb_elem_tot*(nb_som_elem*(nb_som_elem)))+extra_nb_coeff;
  la_matrice.dimensionner(nb_som*nb_comp,nfin*nb_comp,0);
  IntTab Indice(nb_coeff,2);//p(1,2);
  int tot=0;
  for (int elem=0; elem<nb_elem_tot; elem++)
    {
      for (int i1=0; i1<nb_som_elem; i1++)
        {
          int glob=elems(elem,i1);
          for (int i2=0; i2<nb_som_elem; i2++)
            {
              int glob2=elems(elem,i2);
              if (glob<nb_som)
                {
                  Indice(tot,0)=glob;
                  Indice(tot,1)=glob2;
                  tot++;
                }
            }
        }
    }
  // Specific BC that will require extra terms in the matrix:
  for(int i=0; i<le_dom_cl.nb_cond_lim(); i++)
    {
      const Cond_lim_base& ccl = le_dom_cl.les_conditions_limites(i).valeur();
      if(sub_type(Echange_interne_global_impose, ccl))
        {
          const Echange_interne_global_impose& cl = ref_cast(Echange_interne_global_impose, ccl);
          bool is_parfait = sub_type(Echange_interne_global_parfait, ccl);
          const Champ_front_calc_interne& ch = ref_cast(Champ_front_calc_interne, cl.T_ext());
          const IntVect& mp = ch.face_map();   // indices of (local) *faces*
          const Front_VF& le_bord = ref_cast(Front_VF,cl.frontiere_dis());
          int nfacedeb = le_bord.num_premiere_face();
          int nfacefin = nfacedeb + le_bord.nb_faces();
          std::vector<bool> hit(le_bord.nb_faces());
          std::fill(hit.begin(), hit.end(), false);
          for(int face_recto=nfacedeb; face_recto < nfacefin; face_recto++)
            {
              int face_verso = mp(face_recto-nfacedeb)+nfacedeb;
              int face_plus_2 = 0;
 
              if (is_parfait)
                {
                  // element attached to the face
                  int elem1 = le_dom.face_voisins(face_verso, 0);
                  int elem = (elem1 != -1) ? elem1 : le_dom.face_voisins(face_verso, 1);
                  // other face of the element (i.e. face "after" opp_face)
                  int f1 = le_dom.elem_faces(elem, 0);
                  face_plus_2 = f1 != face_verso ? f1 : le_dom.elem_faces(elem, 1);  // 2 faces per elem max - 1D
                }
 
              // all vertices of recto face are put in correspondance with all vertices of verso face:
              for(int som_idx=0; som_idx < nb_som_face; som_idx++)
                {
                  int som_recto = le_dom.face_sommets(face_recto, som_idx);
                  for(int som_idx2=0; som_idx2 < nb_som_face; som_idx2++)
                    {
                      int som_verso = le_dom.face_sommets(face_verso, som_idx2);
                      // tri-diagonal structure extended through the wall
                      // actually this puts too many coeffs in case of is_parfait, but I can't be bothered ...
                      Indice(tot, 0) = som_recto;
                      Indice(tot, 1) = som_verso;
                      tot++;
 
                      if (is_parfait)   // more complicated: recto face takes (i-1), i and ("next face after opp_face")
                        // verso face will just serve to connect recto and verso
                        {
                          int som_plus_2 = le_dom.face_sommets(face_plus_2, som_idx2);
                          if(!hit[face_recto-nfacedeb])  // first inner face
                            {
                              Indice(tot, 0) = som_recto;
                              Indice(tot, 1) = som_plus_2;
                              tot++;
                            }
                          else                  // opposite inner face
                            {
                              // Nothing to do! no extra tri-band coeff.
                            }
                        }
                    }
                }
              hit[face_recto-nfacedeb] = true;
              hit[face_verso-nfacedeb] = true;
            }
        }
    }
 
  Indice.resize(tot,2);
 
  tableau_trier_retirer_doublons(Indice);
  if (nb_comp==1)
    la_matrice.dimensionner(Indice);
  else
    {
      int ntot=Indice.dimension(0);
      // create a nb_comp*nb_comp block for each original coefficient
      IntTab Indice_v(ntot*nb_comp*nb_comp,2);
      int newcas=0;
      for (int cas=0; cas<ntot; cas++)
        {
          for (int i=0; i<nb_comp; i++)
            for (int j=0; j<nb_comp; j++)
              {
                int newi=Indice(cas,0)*nb_comp+i;
                int newj=Indice(cas,1)*nb_comp+j;
                Indice_v(newcas,0)=newi;
                Indice_v(newcas,1)=newj;
                newcas++;
 
              }
        }
      tableau_trier_retirer_doublons(Indice_v);
      la_matrice.dimensionner(Indice_v);
      //Cerr<<Indice<<" nnnnnnnnn "<<Indice_v<<finl;
      // abort();
    }
 
  if (matrice_stocke)
    {
      ref_cast_non_const(Matrice_Morse,matrice_sto_)=la_matrice;
    }
  //
  if (0)
    {
      Cerr<<la_matrice.get_tab1()<<finl;
      Cerr<<la_matrice.get_tab2()<<finl;
      la_matrice.get_set_coeff() = 1;
      Cerr<<"ici "<<tot <<" "<< nb_coeff<<finl;
      Cerr<<"Indice"<<Indice<<finl;
      la_matrice.imprimer_formatte(Cerr);
 
      Process::exit();
    }
 
}
 
/*! @brief Modification of the matrix coefficients and the right-hand side for Dirichlet conditions.
 *
 */
 
void Op_EF_base::modifier_pour_Cl(const Domaine_EF& le_dom,
                                  const Domaine_Cl_EF& le_dom_cl,
                                  Matrice_Morse& la_matrice, DoubleTab& secmem) const
{
  // Sizing the matrix that will receive coefficients from
  // convection and diffusion in the face case.
  // This matrix has a Morse matrix structure.
  // We start by computing the sizes of arrays tab1 and tab2.
  const Conds_lim& les_cl = le_dom_cl.les_conditions_limites();
 
  const auto& tab1 = la_matrice.get_tab1();
  auto& coeff = la_matrice.get_set_coeff();
 
  // determining the size of the unknown field.
  // Cerr << "in modifier cl " << finl;
  const DoubleTab& champ_inconnue = le_dom_cl.equation().inconnue().valeurs();
  int nb_comp = champ_inconnue.line_size();
 
  const IntTab& faces_sommets=le_dom.face_sommets();
  int nb_som_face=faces_sommets.dimension(1);
 
  int nb_som=secmem.dimension(0);
  for (const auto& itr : les_cl)
    {
      const Cond_lim& la_cl = itr;
      const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
      int nfaces = le_bord.nb_faces_tot();
      if (sub_type(Dirichlet,la_cl.valeur()))
        {
          const Dirichlet& la_cl_Dirichlet = ref_cast(Dirichlet,la_cl.valeur());
          for (int ind_face=0; ind_face < nfaces; ind_face++)
            for (int s=0; s<nb_som_face; s++)
              {
                int face=le_bord.num_face(ind_face);
                int som=faces_sommets(face,s);
                if (som<nb_som)
                  for (int comp=0; comp<nb_comp; comp++)
                    {
                      int m=som*nb_comp+comp;
                      auto idiag = tab1[som*nb_comp+comp]-1;
                      auto nbvois = tab1[som*nb_comp+1+comp] - tab1[som*nb_comp+comp];
                      for (auto k=0; k < nbvois; k++)
                        coeff[idiag+k]=0;
                      la_matrice.coef(m,m)=1;
                      assert(la_matrice.coef(som*nb_comp+comp,som*nb_comp+comp)==1);
                      // for neighbors
 
                      // for the right-hand side
                      secmem(som,comp)= la_cl_Dirichlet.val_imp(ind_face,comp);
                      secmem(som,comp)= champ_inconnue(som,comp);
                    }
              }
        }
      if (sub_type(Dirichlet_homogene,la_cl.valeur()))
        {
          const Dirichlet_homogene& la_cl_Dirichlet = ref_cast(Dirichlet_homogene,la_cl.valeur());
          for (int ind_face=0; ind_face < nfaces; ind_face++)
            for (int s=0; s<nb_som_face; s++)
              {
                int face=le_bord.num_face(ind_face);
                int som=faces_sommets(face,s);
 
                for (int comp=0; comp<nb_comp; comp++)
                  {
                    int m=som*nb_comp+comp;
                    auto idiag = tab1[som*nb_comp+comp]-1;
                    auto nbvois = tab1[som*nb_comp+1+comp] - tab1[som*nb_comp+comp];
                    for (auto k=0; k < nbvois; k++)
                      coeff[idiag+k]=0;
                    la_matrice.coef(m,m)=1;
                    assert(la_matrice.coef(som*nb_comp+comp,som*nb_comp+comp)==1);
                    // for neighbors
 
                    // for the right-hand side
                    secmem(som,comp)= 0;
                    if (la_cl_Dirichlet.val_imp(ind_face,comp)!=0)
                      {
                        Cerr<<"not possible, is it ??? "<<finl;
                        Process::exit();
                      }
                  }
              }
        }
 
      if (sub_type(Symetrie,la_cl.valeur()))
        {
          // do nothing here, handled below
        }
    }
 
  // Modify for symmetry
 
  if (le_dom_cl.equation().inconnue().nature_du_champ()==vectoriel)
    {
      le_dom_cl.imposer_symetrie_matrice_secmem(la_matrice,secmem);
    }
 
  /*
    la_matrice.imprimer(Cerr);
    la_matrice.imprimer_formatte(Cerr);
    Cerr<<"scd "<<secmem<<finl;
    */
  //exit();
  // Account for periodic boundary conditions
  //  Cerr << "fin de modifier" << finl;
}
 
/*! @brief Multiplies the boundary flux by rho cp or rho if necessary.
 *
 */
void Op_EF_base::modifier_flux( const Operateur_base& op) const
{
  controle_modifier_flux_=1;
  const Domaine_EF& le_dom_EF=ref_cast(Domaine_EF,op.equation().domaine_dis());
  DoubleTab& flux_bords_=op.flux_bords();
  int nb_compo=flux_bords_.dimension(1);
  const Probleme_base& pb=op.equation().probleme();
 
  // Multiply by rho for incompressible Navier-Stokes
  Nom nom_eqn=op.equation().que_suis_je();
  if (nom_eqn.debute_par("Navier_Stokes") && pb.milieu().que_suis_je()=="Fluide_Incompressible")
    {
      const Champ_base& rho = op.equation().milieu().masse_volumique();
      if (sub_type(Champ_Uniforme,rho))
        {
          double coef = rho.valeurs()(0,0);
          int nb_faces_bord=le_dom_EF.nb_faces_bord();
          for (int face=0; face<nb_faces_bord; face++)
            for(int k=0; k<nb_compo; k++)
              flux_bords_(face,k) *= coef;
        }
    }
 
}
 
/*! @brief Prints the boundary fluxes of an EF operator at faces (i.e., diffusion, convection).
 *
 */
int Op_EF_base::impr(Sortie& os, const Operateur_base& op) const
{
  const Domaine_EF& le_dom_EF=ref_cast(Domaine_EF,op.equation().domaine_dis());
  const DoubleTab& flux_bords_=op.flux_bords();
  if (flux_bords_.nb_dim()!=2)
    {
      Cout << "Flux printing is not implemented for the operator " << op.que_suis_je() << finl;
      return 1;
    }
  if (controle_modifier_flux_==0)
    {
      Cerr<<op.que_suis_je()<<" appelle  Op_EF_base::impr sans avoir appeler  Op_EF_base::modifier_flux, on arrete tout "<<finl;
      Process::exit();
    }
 
  int nb_compo=flux_bords_.dimension(1);
  const Probleme_base& pb=op.equation().probleme();
  const Schema_Temps_base& sch=pb.schema_temps();
  // Print moments only if requested and if we are handling the velocity diffusion operator
  int impr_mom=0;
  if (le_dom_EF.domaine().moments_a_imprimer() && sub_type(Operateur_Diff_base,op) && op.equation().inconnue().le_nom()=="vitesse")
    impr_mom=1;
 
  const int impr_sum=(le_dom_EF.domaine().bords_a_imprimer_sum().est_vide() ? 0:1);
  const int impr_bord=(le_dom_EF.domaine().bords_a_imprimer().est_vide() ? 0:1);
  int flag=0;
  if (Process::je_suis_maitre()) flag=1;
  op.ouvrir_fichier( Flux,"",flag);
  op.ouvrir_fichier(Flux_moment,"moment",impr_mom&&flag);
  op.ouvrir_fichier(Flux_sum,"sum",impr_sum&&flag);
  EcrFicPartage Flux_face;
  op.ouvrir_fichier_partage(Flux_face,"",impr_bord);
 
 
  // Print the current time
  double temps=sch.temps_courant();
  if(Process::je_suis_maitre())
    {
      Flux.add_col(temps);
      if (impr_mom) Flux_moment.add_col(temps);
      if (impr_sum) Flux_sum.add_col(temps);
    }
 
  // Compute moments
  DoubleTab xgr;
  if (impr_mom) xgr = le_dom_EF.calculer_xgr();
  DoubleVect moment(nb_compo);
  moment=0;
 
  // Iterate over boundaries to sum fluxes per boundary into the flux_bord array
  DoubleVect flux_bord(nb_compo);
  DoubleVect bilan(nb_compo);
  bilan = 0;
  for (int num_cl=0; num_cl<le_dom_EF.nb_front_Cl(); num_cl++)
    {
      flux_bord=0;
      const Frontiere_dis_base& la_fr = op.equation().domaine_Cl_dis().les_conditions_limites(num_cl)->frontiere_dis();
      const Cond_lim& la_cl = op.equation().domaine_Cl_dis().les_conditions_limites(num_cl);
      const Front_VF& frontiere_dis = ref_cast(Front_VF,la_cl->frontiere_dis());
      int ndeb = frontiere_dis.num_premiere_face();
      int nfin = ndeb + frontiere_dis.nb_faces();
      for (int face=ndeb; face<nfin; face++)
        {
          for(int k=0; k<nb_compo; k++)
            flux_bord(k)+=flux_bords_(face, k);
 
          if (impr_mom)
            {
              // Compute the moment exerted by the fluid on the boundary (OM/\F)
              if (Objet_U::dimension==2)
                moment(0)+=flux_bords_(face,1)*xgr(face,0)-flux_bords_(face,0)*xgr(face,1);
              else
                {
                  moment(0)+=flux_bords_(face,2)*xgr(face,1)-flux_bords_(face,1)*xgr(face,2);
                  moment(1)+=flux_bords_(face,0)*xgr(face,2)-flux_bords_(face,2)*xgr(face,0);
                  moment(2)+=flux_bords_(face,1)*xgr(face,0)-flux_bords_(face,0)*xgr(face,1);
                }
            }
        }
 
      // Sum contributions from each processor
      for(int k=0; k<nb_compo; k++)
        {
          flux_bord(k)=Process::mp_sum(flux_bord(k));
          if (impr_mom) moment(k)=Process::mp_sum(moment(k));
        }
 
      // Write to output files
      if (Process::je_suis_maitre())
        {
          for(int k=0; k<nb_compo; k++)
            {
              Flux.add_col(flux_bord(k));
              if (impr_mom) Flux_moment.add_col(moment(k));
              if (le_dom_EF.domaine().bords_a_imprimer_sum().contient(la_fr.le_nom())) Flux_sum.add_col(flux_bord(k));
 
              // Sum fluxes from all boundaries into the bilan array
              bilan(k)+=flux_bord(k);
            }
        }
    }
 
  // Print the balance values and move to next line
  if(Process::je_suis_maitre())
    {
      for(int k=0; k<nb_compo; k++)
        Flux.add_col(bilan(k));
 
      Flux << finl;
      if (impr_mom) Flux_moment << finl;
      if (impr_sum) Flux_sum << finl;
    }
 
  // Print per-face values if requested
  for (int num_cl=0; num_cl<le_dom_EF.nb_front_Cl(); num_cl++)
    {
      const Frontiere_dis_base& la_fr = op.equation().domaine_Cl_dis().les_conditions_limites(num_cl)->frontiere_dis();
      const Cond_lim& la_cl = op.equation().domaine_Cl_dis().les_conditions_limites(num_cl);
      const Front_VF& frontiere_dis = ref_cast(Front_VF,la_cl->frontiere_dis());
      int ndeb = frontiere_dis.num_premiere_face();
      int nfin = ndeb + frontiere_dis.nb_faces();
      // Print per-face values
      if (le_dom_EF.domaine().bords_a_imprimer().contient(la_fr.le_nom()))
        {
          Flux_face << "# Flux par face sur " << la_fr.le_nom() << " au temps " << temps << " : " << finl;
          const DoubleTab& xv=le_dom_EF.xv();
          for (int face=ndeb; face<nfin; face++)
            {
              if (Objet_U::dimension==2)
                Flux_face << "# Face a x= " << xv(face,0) << " y= " << xv(face,1) << " flux=" ;
              else if (Objet_U::dimension==3)
                Flux_face << "# Face a x= " << xv(face,0) << " y= " << xv(face,1) << " z= " << xv(face,2) << " flux=" ;
              for(int k=0; k<nb_compo; k++)
                Flux_face << " " << flux_bords_(face, k);
              Flux_face << finl;
            }
          Flux_face.syncfile();
        }
    }
 
  return 1;
}
 
int Op_EF_base::elem_contribue(const int elem) const
{
  if (marqueur_elem_.size_array()==0)
    return 1;
  else if (marqueur_elem_[elem]==1)
    return 0;
  return 1;
}
 
void Op_EF_base::marque_elem(const Equation_base& eqn)
{
  if (eqn.has_champ("marqueur_loi_de_paroi"))
    {
      const DoubleTab& marq = eqn.get_champ("marqueur_loi_de_paroi").valeurs();
      int ntot = marq.dimension_tot(0);
      marqueur_elem_.resize_array(ntot);
      for (int n = 0; n < ntot; n++)
        if (marq(n) > 0)
          marqueur_elem_[n] = 1;
    }
}