Champ_Post_Operateur_Eqn.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_Post_Operateur_Eqn.h>
#include <Probleme_base.h>
#include <Postraitement.h>
#include <Domaine_VF.h>
#include <Nom.h>
#include <Champ_Generique_refChamp.h>
#include <Discretisation_base.h>
#include <Equation_base.h>
#include <Operateur_base.h>
#include <Param.h>
#include <MD_Vector_composite.h>
 
Implemente_instanciable_sans_constructeur(Champ_Post_Operateur_Eqn,"Operateur_Eqn|Champ_Post_Operateur_Eqn",Champ_Generique_Operateur_base);
// XD champ_post_operateur_eqn champ_post_de_champs_post operateur_eqn BRACE Post-process equation operators/sources
 
Sortie& Champ_Post_Operateur_Eqn::printOn(Sortie& s ) const
{
  return s << que_suis_je() << " " << le_nom();
}
Champ_Post_Operateur_Eqn::Champ_Post_Operateur_Eqn()
{
}
 
void Champ_Post_Operateur_Eqn::set_param(Param& param) const
{
  Champ_Generique_Operateur_base::set_param(param);
  param.ajouter("numero_source",&numero_source_); // XD_ADD_P entier
  // XD_CONT the source to be post-processed (its number). If you have only one source term, numero_source will
  // XD_CONT correspond to 0 if you want to post-process that unique source
  param.ajouter("numero_op",&numero_op_); // XD_ADD_P entier
  // XD_CONT numero_op will be 0 (diffusive operator) or 1 (convective operator) or 2 (gradient operator) or 3
  // XD_CONT (divergence operator).
  param.ajouter("numero_masse",&numero_masse_); // XD_ADD_P entier
  // XD_CONT numero_masse will be 0 for the mass equation operator in Pb_multiphase.
  param.ajouter_flag("sans_solveur_masse",&sans_solveur_masse_); // XD_ADD_P rien
  // XD_CONT not_set
  param.ajouter("compo",&compo_); // XD_ADD_P entier
  // XD_CONT If you want to post-process only one component of a vector field, you can specify the number of the
  // XD_CONT component after compo keyword. By default, it is set to -1 which means that all the components will be
  // XD_CONT post-processed. This feature is not available in VDF disretization.
}
 
Entree& Champ_Post_Operateur_Eqn::readOn(Entree& s )
{
  Champ_Generique_Operateur_base::readOn(s);
 
  return s ;
}
 
void Champ_Post_Operateur_Eqn::verification_cas_compo() const
{
  // Applying compo to a vector
  if ((!ref_eq_->inconnue().is_vectorial()) && (compo_ != -1 ))
    {
      Cerr<<"Error in Champ_Post_Operateur_Eqn::verification_cas_compo()"<<finl;
      Cerr<<"It isn't possible to get a component from a non vectoriel field " <<finl;
      Process::exit();
    }
 
  // Checking compo
  const int nb_compo= ref_eq_->inconnue().nb_comp();
  if ((compo_<-1)||(compo_>nb_compo-1))
    {
      Cerr<<"Error in Champ_Post_Operateur_Eqn::verification_cas_compo()"<<finl;
      Cerr<<"compo="<<compo_<<" is not allowed."<<" You must give a value between 0 and "<< nb_compo-1<<finl;
      Process::exit();
    }
  // Check that we are not in VDF
  if (ref_eq_->discretisation().is_vdf() && (compo_ != -1 ))
    {
      Cerr<<"Error in Champ_Post_Operateur_Eqn::verification_cas_compo()"<<finl;
      Cerr<<"The option compo is not available in case of VDF discretization"<<finl;
      Process::exit();
    }
}
 
void Champ_Post_Operateur_Eqn::completer(const Postraitement_base& post)
{
 
  Champ_Gen_de_Champs_Gen::completer(post);
 
  const Probleme_base& Pb = ref_cast(Postraitement,post).probleme();
  int numero_eq_=-1;
  if (sub_type(Champ_Generique_refChamp,get_source(0)))
    {
 
      OWN_PTR(Champ_base) espace_stockage;
      const Champ_base& mon_champ = get_source(0).get_champ(espace_stockage);
      if (sub_type(Champ_Inc_base,mon_champ))
        {
          const Champ_Inc_base& mon_champ_inc = ref_cast(Champ_Inc_base,mon_champ);
 
          {
 
            //Retrieving the equation before it is associated with the Champ_Inc
            //Iterating over the problem equations and identifying the one corresponding to the unknown field
 
            int nb_eq = Pb.nombre_d_equations();
            int i=0;
 
            while (i<nb_eq)
              {
                const Equation_base& eq_test = Pb.equation(i);
                if ((eq_test.inconnue().le_nom() == mon_champ_inc.le_nom()))
                  {
                    numero_eq_=i;
                    break;
                  }
                i++;
              }
          }
        }
    }
  if (numero_eq_==-1)
    {
      Cerr<<"Champ_Post_Operateur_Eqn can be applied only to equation unknown."<<finl;
      Process::exit();
    }
 
  ref_eq_=Pb.equation(numero_eq_);
 
  int ok=0;
  const Equation_base& eqn=ref_eq_.valeur();
  const MD_Vector& mdf = eqn.inconnue().valeurs().get_md_vector(),
                   md = sub_type(MD_Vector_composite, mdf.valeur()) ? ref_cast(MD_Vector_composite, mdf.valeur()).get_desc_part(0) : mdf;
  const Domaine_VF& zvf= ref_cast( Domaine_VF,ref_eq_->domaine_dis());
  if (md== zvf.face_sommets().get_md_vector())
    {
      localisation_inco_=Entity::FACE;
      ok=1;
    }
  if (md== zvf.domaine().les_elems().get_md_vector())
    {
      localisation_inco_=Entity::ELEMENT;
      ok=1;
    }
  if (md == zvf.domaine().les_sommets().get_md_vector())
    {
      ok=1;
      localisation_inco_=Entity::NODE;
    }
  if (ok==0)
    {
      Cerr<<"Error in "<<que_suis_je()<<" unknown localisation"<<finl;
      Process::exit();
    }
  verification_cas_compo();
}
 
const Champ_base& Champ_Post_Operateur_Eqn::get_champ_without_evaluation(OWN_PTR(Champ_base)& espace_stockage) const
{
  espace_stockage = ref_eq_->inconnue();
  return espace_stockage;
}
 
const Champ_base& Champ_Post_Operateur_Eqn::get_champ_compo_without_evaluation(OWN_PTR(Champ_base)& espace_stockage) const
{
 
  OWN_PTR(Champ_Fonc_base)  espace_stockage_fonc;
  //Champ source_espace_stockage;
 
  double temps=0.;
  Nom directive;
 
  switch (localisation_inco_)
    {
    case Entity::ELEMENT:
      directive="CHAMP_ELEM";
      break;
    case Entity::NODE:
      directive="CHAMP_SOMMETS";
      break;
    case Entity::FACE:
      directive="CHAMP_FACE";
      break;
    default:
      Cerr<<"error in Champ_Post_Operateur_Eqn::get_champ"<<finl;
      Process::exit();
 
    }
  int nb_comp = 1;
  ref_eq_->discretisation().discretiser_champ(directive,ref_eq_->domaine_dis(),"oooo","unit", nb_comp,temps,espace_stockage_fonc);
  espace_stockage = espace_stockage_fonc;
  espace_stockage->fixer_nature_du_champ(scalaire);
 
  return espace_stockage;
}
 
const Champ_base& Champ_Post_Operateur_Eqn::get_champ(OWN_PTR(Champ_base)& espace_stockage) const
{
  // Start by building the complete vector field
  OWN_PTR(Champ_base) espace_stockage_complet;
  espace_stockage_complet = get_champ_without_evaluation(espace_stockage_complet);
  DoubleTab& es = (espace_stockage_complet->valeurs());
 
  //if (ref_eq_->schema_temps().temps_courant()!=0)
  {
    if (numero_op_!=-1)
      {
        // some compute methods are wrong!!!! they should all be recoded as res=0 ajouter();
        es=0;
        Operateur().ajouter(ref_eq_->operateur(numero_op_).mon_inconnue().valeurs(),es);
      }
    else if (numero_source_!=-1)
      ref_eq_->sources()(numero_source_).calculer(es);
    else if ((numero_masse_!=-1) && ref_eq_->has_interface_blocs())
      es=0, ref_eq_->schema_temps().ajouter_blocs({},es,ref_eq_.valeur());
    if (!sans_solveur_masse_)
      ref_eq_->solv_masse().appliquer_impl(es); //We divide by the volume
    // Hack: because Masse_PolyMAC_CDO_Face::appliquer_impl does not divide by the volume (mass matrix)....
    if (ref_eq_->solv_masse().que_suis_je()=="Masse_PolyMAC_CDO_Face")
      {
        //Cerr << "Volumic source terms on faces with PolyMAC_CDO can't be post-processed yet." << finl;
        Cerr << "Warning, source terms on faces with PolyMAC_CDO are post-processed as S*dV not as volumic source terms S." << finl;
        Cerr << "Cause Masse_PolyMAC_CDO_Face::appliquer_impl do not divide per volume." << finl;
        //Process::exit();
      }
  }
  es.echange_espace_virtuel();
  // then I retrieve only the component if it was requested
  if (compo_>-1)
    {
      // prepare the storage space for a single component
      OWN_PTR(Champ_base) espace_stockage_compo;
      espace_stockage_compo = get_champ_compo_without_evaluation(espace_stockage_compo);
      DoubleTab& es_compo = (espace_stockage_compo->valeurs());
      int nb_pos = es.dimension(0);
      for (int i=0; i<nb_pos; i++)
        {
          es_compo(i) = es(i,compo_);
        }
      espace_stockage = espace_stockage_compo;
    }
  else
    espace_stockage = espace_stockage_complet;
 
  return espace_stockage;
}
 
const Noms Champ_Post_Operateur_Eqn::get_property(const Motcle& query) const
{
//Component creation should be done dynamically (Eqn_...)
 
  Motcles motcles(2);
  motcles[0] = "composantes";
  motcles[1] = "unites";
 
  int rang = motcles.search(query);
 
  switch(rang)
    {
    case 0:
      {
        if (compo_==-1)
          {
            int nb_comp= ref_eq_->inconnue().nb_comp();
            Noms compo(nb_comp);
            for (int i=0; i<nb_comp; i++)
              {
                Nom nume(i);
                compo[i] = nom_post_+nume;
              }
            return compo;
          }
        else
          {
            Noms compo(1);
            compo[0] = nom_post_;
            return compo;
          }
      }
 
    case 1:
      {
        if (compo_==-1)
          {
            int nb_comp= ref_eq_->inconnue().nb_comp();
            Noms unites(nb_comp);
            //Noms source_unites = get_source(0).get_property("unites");
            for (int i=0; i<nb_comp; i++)
              {
                unites[i] = "unit";
              }
            return unites;
          }
        else
          {
            // Using a vector because the method returns Noms
            Noms unites(1);
            unites[0] = "unit";
            return unites;
          }
      }
 
    }
  return Champ_Gen_de_Champs_Gen::get_property(query);
}
 
Entity Champ_Post_Operateur_Eqn::get_localisation(const int index) const
{
 
  return localisation_inco_;
 
}
//Name the field as a source by default
//"Eqn_" + source_field_name
void Champ_Post_Operateur_Eqn::nommer_source()
{
  if (nom_post_=="??")
    {
      Nom nom_post_source;
      //const Noms nom = get_source(0).get_property("nom");
      nom_post_source = "Eqn_s";
      nom_post_source += Nom(numero_source_);
      nom_post_source += "_o" ;
      nom_post_source += Nom(numero_op_);
      nom_post_source += "_m" ;
      nom_post_source += Nom(numero_masse_);
      if (compo_!=-1)
        {
          Nom nume(compo_);
          nom_post_source += nume;
        }
      nommer(nom_post_source);
    }
 
}
const Operateur_base& Champ_Post_Operateur_Eqn::Operateur() const
{
  return ref_eq_->operateur(numero_op_).l_op_base();
}
 
Operateur_base& Champ_Post_Operateur_Eqn::Operateur()
{
  return ref_eq_->operateur(numero_op_).l_op_base();
}