Op_Correlation.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_Generique_refChamp.h>
#include <Discretisation_base.h>
#include <Schema_Temps_base.h>
#include <Op_Correlation.h>
 
Implemente_instanciable(Op_Correlation, "Op_Correlation", Operateur_Statistique_tps_base);
 
Sortie& Op_Correlation::printOn(Sortie& s) const
{
  return s << que_suis_je() << " " << le_nom();
}
 
Entree& Op_Correlation::readOn(Entree& s)
{
  return s;
}
 
void Op_Correlation::completer(const Probleme_base& Pb, const Nom& prefix)
{
  const OBS_PTR(Champ_Generique_base)& a = la_moyenne_a().integrale().le_champ();
  const OBS_PTR(Champ_Generique_base)& b = la_moyenne_b().integrale().le_champ();
  integrale_tps_a_ = la_moyenne_a().integrale();
  integrale_tps_b_ = la_moyenne_b().integrale();
 
  // Retrieve the discretized domain
  const Domaine_dis_base& domaine = a->get_ref_domaine_dis_base();
 
  // Size the integrale_champ field to the same size as mon_champ
  const DoubleTab& tab_a = valeurs_a();
  const DoubleTab& tab_b = valeurs_b();
 
  OWN_PTR(Champ_base) espace_stockage_source_a, espace_stockage_source_b;
  const Champ_base& source_a = a->get_champ(espace_stockage_source_a);
  const Champ_base& source_b = b->get_champ(espace_stockage_source_b);
 
  const int nb_compo_a = source_a.nb_comp();
  const int nb_compo_b = source_b.nb_comp();
 
  // Use nom_cible only if both a and b are Champ_Generique_refChamp,
  // so that the variable names can be reproduced in the lml output.
  Noms noms_a, noms_b;
 
  if (!(sub_type(Champ_Generique_refChamp,a.valeur()) && (sub_type(Champ_Generique_refChamp, b.valeur()))))
    {
      noms_a = a->get_property("nom");
      noms_b = b->get_property("nom");
    }
  else
    {
      noms_a = a->get_property("nom_cible");
      noms_b = b->get_property("nom_cible");
    }
 
  const Nom nom_a = noms_a[0];
  const Nom nom_b = noms_b[0];
 
  // Special case when a and b are Champ_Generique_refChamp,
  // so that the variable names can be reproduced in the lml output.
  Noms compo_a, compo_b;
 
  if (!(sub_type(Champ_Generique_refChamp,a.valeur()) && (sub_type(Champ_Generique_refChamp, b.valeur()))))
    {
      compo_a = a->get_property("composantes");
      compo_b = b->get_property("composantes");
    }
  else
    {
      compo_a = source_a.noms_compo();
      compo_b = source_b.noms_compo();
    }
 
  const Noms unites_a = a->get_property("unites");
  const Noms unites_b = b->get_property("unites");
 
  int nb_valeurs_nodales = -1;
 
  // Note B.M.: integrale_tps_ab_ is already typed in Correlation::associer,
  // but fixer_nb_valeurs_nodales has not yet been called.
 
  if ((tab_a.get_md_vector() != tab_b.get_md_vector())
      || (tab_a.dimension_tot(0) != tab_b.dimension_tot(0))
      || (nb_compo_a > 1 && (tab_a.nb_dim() == 1 || tab_a.dimension(1) != nb_compo_a))
      || (nb_compo_b > 1 && (tab_b.nb_dim() == 1 || tab_b.dimension(1) != nb_compo_b)))
    {
      // Different support for fields a and b or for their components:
      // bring everything to elements by creating a Champ_Fonc_P0_VDF or VEF.
      integrale_tps_ab_.support_different() = 1;
      Nom Champ_Fonc_P0 = "Champ_Fonc_P0_";
      Champ_Fonc_P0 += Pb.discretisation().que_suis_je().substr_old(1, 3);
      integrale_tps_ab_.typer_champ(Champ_Fonc_P0);
      nb_valeurs_nodales = domaine.nb_elem();
    }
  else
    {
      // Same support for fields a and b
      integrale_tps_ab_.support_different() = 0;
      if (tab_a.size_reelle_ok())
        nb_valeurs_nodales = tab_a.dimension(0);
      else
        nb_valeurs_nodales = -1; // B.M.: quoi donner a fixer_nb_valeurs_nodales ???
    }
  int nb_comp = nb_compo_a * nb_compo_b;
  DoubleTab& val = valeurs_ab();
  Noms noms_composantes;
  if (nb_compo_a == 1)
    if (nb_compo_b == 1)
      {
        // A and B are scalar fields
      }
    else
      {
        // A is a scalar field, B is a vector field
        noms_composantes.dimensionner(nb_comp);
        Nom debut("Correlation_");
        debut += nom_a + "_";
        for (int i = 0; i < nb_comp; i++)
          noms_composantes[i] = debut + compo_b[i];
      }
  else if (nb_compo_b == 1)
    {
      // A is a vector field, B is a scalar field
      noms_composantes.dimensionner(nb_comp);
      Nom debut("Correlation_");
      debut += nom_b + "_";
      for (int i = 0; i < nb_comp; i++)
        noms_composantes[i] = debut + compo_a[i];
    }
  else
    {
      // A and B are vector fields
      noms_composantes.dimensionner(nb_comp);
      Nom debut("Correlation_");
      for (int i = 0; i < nb_compo_a; i++)
        for (int j = 0; j < nb_compo_b; j++)
          noms_composantes[i * nb_compo_b + j] = debut + compo_a[i] + "_" + compo_b[j];
    }
 
  if (nb_comp > 1)
    val.resize(0, nb_comp);
 
  Nom nom_pour_post = "Correlation_";
  nom_pour_post += nom_a + "_" + nom_b;
  integrale_tps_ab_.le_champ_calcule().nommer(nom_pour_post);
 
  Nom pdi_name = prefix + nom_pour_post;
  integrale_tps_ab_.le_champ_calcule().set_pdi_name(pdi_name);
 
  integrale_tps_ab_.le_champ_calcule().associer_domaine_dis_base(domaine);
  integrale_tps_ab_.le_champ_calcule().fixer_nb_comp(nb_comp);
 
  // BM: le parametre de fixer_nb_valeurs_nodales est inutile,
  //  The field can determine on its own how many DOFs it has, but
  //  there is a check on the parameter: the value must be correct.
  integrale_tps_ab_.le_champ_calcule().fixer_nb_valeurs_nodales(nb_valeurs_nodales);
 
  // Set the units
  Nom unite(unites_a[0]);
  unite += ".";
  unite += unites_b[0];
 
  if (nb_comp == 1)
    integrale_tps_ab_.le_champ_calcule().fixer_unite(unite);
  else
    {
      Noms unites(nb_comp);
      for (int i = 0; i < unites.size(); i++)
        unites[i] = unite;
      integrale_tps_ab_.le_champ_calcule().fixer_noms_compo(noms_composantes);
      integrale_tps_ab_.le_champ_calcule().fixer_unites(unites);
    }
  integrale_tps_ab_.le_champ_calcule().changer_temps(Pb.schema_temps().temps_courant());
}
 
DoubleTab Op_Correlation::calculer_valeurs() const
{
  Integrale_tps_produit_champs correlation(integrale_tps_ab_);
  const double dt_ab = dt_integration_ab();
  if (dt_ab > 0)
    {
      // Compute Mean(a'b') = Mean(ab) - Mean(a)*Mean(b)
      correlation.le_champ_calcule().valeurs() /= dt_ab;
      const double dt_a = dt_integration_a();
      const double dt_b = dt_integration_b();
      assert(est_egal(dt_a, dt_ab));
      assert(est_egal(dt_b, dt_ab));
      correlation.ajoute_produit_tensoriel(-1 / (dt_a * dt_b), integrale_tps_a_->le_champ_calcule(), integrale_tps_b_->le_champ_calcule());
    }
  return correlation.le_champ_calcule().valeurs();
}
 
int Op_Correlation::completer_post_statistiques(const Domaine& dom, const int is_axi, Format_Post_base& format)
{
  return 1;
}