VEF_discretisation.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 <Taux_cisaillement_P0_VEF.h>
#include <Champ_Fonc_Tabule_P0_VEF.h>
#include <Correlation_Vec_Sca_VEF.h>
#include <Rotationnel_Champ_P1NC.h>
#include <Rotationnel_Champ_Q1NC.h>
#include <Critere_Q_Champ_P1NC.h>
#include <T_paroi_Champ_P1NC.h>
#include <VEF_discretisation.h>
#include <Y_plus_Champ_P1NC.h>
#include <Navier_Stokes_std.h>
#include <grad_T_Champ_P1NC.h>
#include <grad_U_Champ_P1NC.h>
#include <h_conv_Champ_P1NC.h>
#include <Champ_Ostwald_VEF.h>
#include <Champ_Fonc_Tabule.h>
#include <Schema_Temps_base.h>
#include <Champ_Fonc_Q1NC.h>
#include <Fluide_Ostwald.h>
#include <Champ_Uniforme.h>
#include <Domaine_VEF.h>
#include <Quadri_VEF.h>
#include <Tetra_VEF.h>
#include <Hexa_VEF.h>
#include <Tri_VEF.h>
 
Implemente_instanciable(VEF_discretisation, "VEFPreP1B|VEF", Discret_Thyd);
// XD vef discretisation_base vefprep1b INHERITS_BRACE Finite element volume discretization (P1NC/P1-bubble element).
// XD_CONT Since the 1.5.5 version, several new discretizations are available thanks to the optional keyword Read. By
// XD_CONT default, the VEFPreP1B keyword is equivalent to the former VEFPreP1B formulation (v1.5.4 and sooner). P0P1
// XD_CONT (if used with the strong formulation for imposed pressure boundary) is equivalent to VEFPreP1B but the
// XD_CONT convergence is slower. VEFPreP1B dis is equivalent to VEFPreP1B dis Read dis { P0 P1
// XD_CONT Changement_de_base_P1Bulle 1 Cl_pression_sommet_faible 0 }
 
// By default, P0+P1 and changing basis to do P0+P1->P1Bulle
// REMEMBER: we don't necessarily go through the readOn() piece of code if we just write
//     VEF dis
// instead of
//     VEF dis
//     Read dis { ... }
// !!!!
Entree& VEF_discretisation::readOn(Entree& is)
{
  alphaE_ = false;
  alphaS_ = false;
  P1Bulle_ = 0;
  cl_pression_sommet_faible_ = 1;
  return Discret_Thyd::readOn(is);
}
 
void VEF_discretisation::set_param(Param& param) const
{
  Discret_Thyd::set_param(param);
 
  param.ajouter("changement_de_base_P1bulle", &P1Bulle_); // XD_ADD_P entier(into=[0,1])
  // XD_CONT changement_de_base_p1bulle 1 This option may be used to have the P1NC/P0P1 formulation (value set to 0) or
  // XD_CONT the P1NC/P1Bulle formulation (value set to 1, the default).
  param.ajouter_flag("P0", &alphaE_); // XD_ADD_P rien
  // XD_CONT Pressure nodes are added on element centres
  param.ajouter_flag("P1", &alphaS_); // XD_ADD_P rien
  // XD_CONT Pressure nodes are added on vertices
  param.ajouter_flag("Pa", &alphaA_); // XD_ADD_P rien
  // XD_CONT Only available in 3D, pressure nodes are added on bones
  param.ajouter_flag("RT", &alphaRT_); // XD_ADD_P rien
  // XD_CONT For P1NCP1B (in TrioCFD)
  param.ajouter("modif_div_face_dirichlet", &modif_div_face_dirichlet_); // XD_ADD_P entier(into=[0,1])
  // XD_CONT This option (by default 0) is used to extend control volumes for the momentum equation.
  param.ajouter("CL_pression_sommet_faible", &cl_pression_sommet_faible_); // XD_ADD_P entier(into=[0,1])
  // XD_CONT This option is used to specify a strong formulation (value set to 0, the default) or a weak formulation
  // XD_CONT (value set to 1) for an imposed pressure boundary condition. The first formulation converges quicker and is
  // XD_CONT stable in general cases. The second formulation should be used if there are several outlet boundaries with
  // XD_CONT Neumann condition (see Ecoulement_Neumann test case for example).
}
 
void VEF_discretisation::check_param()
{
  // A few checks
  if (dimension != 3 && alphaA_)
    Process::exit("Pa support is only available in 3D.");
  if (!alphaE_ && !alphaS_ && !alphaA_)
    {
      // P0P1Bulle by default with strong pressure BC
      alphaE_ = true;
      alphaS_ = true;
      P1Bulle_ = true;
      cl_pression_sommet_faible_ = 0;
      //Process::exit("You must choose at least one support among P0, P1, Pa.");
    }
  if (alphaA_ && !alphaE_)
    Process::exit("Discretisations Pa or P1+Pa are not yet supported.");
  if (P1Bulle_)
    if (!((alphaE_) && (alphaS_) && (!alphaA_)))
      Process::exit("Option 'changement_de_base_P1bulle' is available only in P0/P1");
 
  if (alphaRT_)
    {
      if (!alphaE_ && !alphaS_ && !alphaA_)
        Process::exit("Choose P0 discretization.");
      Cerr << "Linke's scheme, OK for steady-state pb. The pressure computed is the Bernoulli pressure P =p+|u|^2" << finl;
    }
}
 
Sortie& VEF_discretisation::printOn(Sortie& s) const
{
  s << "lire " << le_nom() << " { " << finl;
  if (alphaE_) s << " P0";
  if (alphaS_) s << " P1";
  if (alphaA_) s << " Pa";
  if (alphaRT_) s << " RT";
  s << " changement_de_base_P1bulle " << P1Bulle_;
  s << " modif_div_face_dirichlet " << modif_div_face_dirichlet_;
  s << " cl_pression_sommet_faible " << cl_pression_sommet_faible_;
  s << " } " << finl;
  return s;
}
 
/*! @brief Discretizes a field for VEFP1B according to a discretization directive.
 *
 * @brief The directive is a Motcle such as "pression",
 *  "divergence_vitesse" (creates a field of type P1_isoP1Bulle).
 *  This method determines the field type to create based on element type
 *  and the discretization directive. It then determines the number of degrees of freedom
 *  and sets all field parameters (type, nb_compo, nb_ddl, nb_pas_dt,
 *  name(s), unit(s), field nature, and assigns a time) and associates the Domaine_dis to the field.
 *  See the code for the correspondence between directives and
 *  the field type created.
 *
 */
void VEF_discretisation::discretiser_champ(const Motcle& directive, const Domaine_dis_base& z, Nature_du_champ nature, const Noms& noms, const Noms& unites, int nb_comp, int nb_pas_dt, double temps,
                                           OWN_PTR(Champ_Inc_base)& champ, const Nom& sous_type) const
{
  Motcles motcles(2);
  motcles[0] = "pression";    // Standard choice for pressure
  motcles[1] = "divergence_vitesse"; // Field type obtained by computing div v
 
  Nom type;
  int default_nb_comp = 0; // Default number of components
  int rang = motcles.search(directive);
  switch(rang)
    {
    case 0:
      type = "Champ_P1_isoP1Bulle";
      default_nb_comp = 1;
      break;
    case 1:
      type = "Champ_P1_isoP1Bulle";
      default_nb_comp = 1;
      break;
    default:
      assert(rang < 0);
      break;
    }
 
  if (directive == DEMANDE_DESCRIPTION)
    Cerr << "VEF_discretisation : " << motcles;
 
  // If the directive is not recognized (or if it is a description request)
  // then call the ancestor:
  if (rang < 0)
    {
      VEF_discretisation::discretiser_champ_(directive, z, nature, noms, unites, nb_comp, nb_pas_dt, temps, champ);
      return;
    }
 
  const int nb_ddl = -1; // the p1b descriptor will provide the actual number
  if (nb_comp < 0)
    nb_comp = default_nb_comp;
  assert(nb_comp > 0);
  creer_champ(champ, z, type, noms[0], unites[0], nb_comp, nb_ddl, nb_pas_dt, temps, directive, que_suis_je());
 
  if (nature == multi_scalaire)
    {
      Cerr << "There is no field of type OWN_PTR(Champ_Inc_base) with P1Bulle discretization" << finl;
      Cerr << "and a multi_scalaire nature." << finl;
      exit();
    }
}
 
void VEF_discretisation::discretiser_champ_(const Motcle& directive, const Domaine_dis_base& z, Nature_du_champ nature, const Noms& noms, const Noms& unites, int nb_comp, int nb_pas_dt, double temps,
                                            OWN_PTR(Champ_Inc_base)& champ, const Nom& sous_type) const
{
  const Domaine_VEF& domaine_vef = ref_cast(Domaine_VEF, z);
 
  Motcles motcles(7);
  motcles[0] = "vitesse";     // Standard choice for velocity
  motcles[1] = "pression";    // Standard choice for pressure
  motcles[2] = "temperature"; // Standard choice for temperature
  motcles[3] = "divergence_vitesse"; // Field type obtained by computing div v
  motcles[4] = "gradient_pression";  // Field type obtained by computing grad P
  motcles[5] = "champ_elem";    // Create a field at elements (type P0)
  motcles[6] = "champ_sommets"; // Create a field at vertices (type P1)
 
  // The velocity field type depends on element type:
  Nom type_champ_vitesse;
  if (sub_type(Tri_VEF, domaine_vef.type_elem()) || sub_type(Tetra_VEF, domaine_vef.type_elem()))
    type_champ_vitesse = "Champ_P1NC";
  else if (sub_type(Quadri_VEF,domaine_vef.type_elem()) || sub_type(Hexa_VEF, domaine_vef.type_elem()))
    type_champ_vitesse = "Champ_Q1NC";
  else
    {
      Cerr << "VEF_discretisation::discretiser_champ :\n The geometric element ";
      Cerr << domaine_vef.type_elem().que_suis_je();
      Cerr << " is not supported." << finl;
      exit();
    }
 
  Nom type;
  int default_nb_comp = 0; // Default number of components
  int rang = motcles.search(directive);
  switch(rang)
    {
    case 0:
      type = type_champ_vitesse;
      default_nb_comp = dimension;
      break;
    case 1:
      type = "Champ_P0_VEF";
      default_nb_comp = 1;
      break;
    case 2:
      type = type_champ_vitesse;
      default_nb_comp = 1;
      break;
    case 3:
      type = "Champ_P0_VEF";
      default_nb_comp = 1;
      break;
    case 4:
      type = type_champ_vitesse;
      default_nb_comp = dimension;
      break;
    case 5:
      type = "Champ_P0_VEF";
      default_nb_comp = 1;
      break;
    case 6:
      type = "Champ_P1_VEF";
      default_nb_comp = 1;
      break;
    default:
      assert(rang < 0);
      break;
    }
 
  if (directive == DEMANDE_DESCRIPTION)
    Cerr << "VEF_discretisation : " << motcles;
 
  // If the directive is not recognized (or if it is a description request)
  // then call the ancestor:
  if (rang < 0)
    {
      Discret_Thyd::discretiser_champ(directive, z, nature, noms, unites, nb_comp, nb_pas_dt, temps, champ);
      return;
    }
 
  // Compute the number of degrees of freedom
  int nb_ddl = 0;
  if (type == "Champ_P0_VEF")
    nb_ddl = z.nb_elem();
  else if (type == type_champ_vitesse)
    nb_ddl = domaine_vef.nb_faces();
  else if (type == "Champ_P1_VEF")
    nb_ddl = domaine_vef.nb_som();
  else
    assert(0);
 
  if (nb_comp < 0)
    nb_comp = default_nb_comp;
  assert(nb_comp > 0);
  creer_champ(champ, z, type, noms[0], unites[0], nb_comp, nb_ddl, nb_pas_dt, temps, directive, que_suis_je());
  if (nature == multi_scalaire)
    {
      champ->fixer_nature_du_champ(nature);
      champ->fixer_unites(unites);
      champ->fixer_noms_compo(noms);
    }
 
}
 
/*! @brief Same as VEF_discretisation::discretiser_champ(..., Champ_Inc).
 *
 */
void VEF_discretisation::discretiser_champ(const Motcle& directive, const Domaine_dis_base& z, Nature_du_champ nature, const Noms& noms, const Noms& unites, int nb_comp, double temps,
                                           OWN_PTR(Champ_Fonc_base)& champ) const
{
  discretiser_champ_fonc_don(directive, z, nature, noms, unites, nb_comp, temps, champ);
}
 
/*! @brief Same as VEF_discretisation::discretiser_champ(..., Champ_Inc).
 *
 */
void VEF_discretisation::discretiser_champ(const Motcle& directive, const Domaine_dis_base& z, Nature_du_champ nature, const Noms& noms, const Noms& unites, int nb_comp, double temps,
                                           OWN_PTR(Champ_Don_base)& champ) const
{
  discretiser_champ_fonc_don(directive, z, nature, noms, unites, nb_comp, temps, champ);
}
 
/*! @brief Same as VEF_discretisation::discretiser_champ(..., Champ_Inc). Common processing for champ_fonc and champ_don.
 *  This method is private (passing an Objet_U is not clean from the outside).
 *
 */
void VEF_discretisation::discretiser_champ_fonc_don(const Motcle& directive, const Domaine_dis_base& z, Nature_du_champ nature, const Noms& noms, const Noms& unites, int nb_comp, double temps,
                                                    Objet_U& champ) const
{
  // Two pointers for easy access to champ_don or champ_fonc depending on the type of the champ object.
  OWN_PTR(Champ_Fonc_base)  * champ_fonc = dynamic_cast<OWN_PTR(Champ_Fonc_base)*>(&champ);
  OWN_PTR(Champ_Don_base) * champ_don = dynamic_cast<OWN_PTR(Champ_Don_base)*>(&champ);
 
  Motcles motcles(2);
  motcles[0] = "pression";    // Standard choice for pressure
  motcles[1] = "divergence_vitesse"; // Field type obtained by computing div v
 
  Nom type;
  int default_nb_comp = 0; // Default number of components
  int rang = motcles.search(directive);
  switch(rang)
    {
    case 0:
      type = "Champ_Fonc_P1_isoP1Bulle";
      default_nb_comp = 1;
      break;
    case 1:
      type = "Champ_Fonc_P1_isoP1Bulle";
      default_nb_comp = 1;
      break;
    default:
      assert(rang < 0);
      break;
    }
 
  if (directive == DEMANDE_DESCRIPTION)
    Cerr << "VEF_discretisation : " << motcles;
 
  // If the directive is not recognized (or if it is a description request)
  // then call the ancestor:
  if (rang < 0)
    {
      if (champ_fonc)
        VEF_discretisation::discretiser_champ_fonc_don_(directive, z, nature, noms, unites, nb_comp, temps, *champ_fonc);
      else
        VEF_discretisation::discretiser_champ_fonc_don_(directive, z, nature, noms, unites, nb_comp, temps, *champ_don);
      return;
    }
 
  int nb_ddl = -1; // the p1b descriptor will provide the actual number
 
  if (nb_comp < 0)
    nb_comp = default_nb_comp;
  assert(nb_comp > 0);
  if (champ_fonc)
    creer_champ(*champ_fonc, z, type, noms[0], unites[0], nb_comp, nb_ddl, temps, directive, que_suis_je());
  else
    creer_champ(*champ_don, z, type, noms[0], unites[0], nb_comp, nb_ddl, temps, directive, que_suis_je());
 
  if (nature == multi_scalaire)
    {
      Cerr << "There is no field of type OWN_PTR(Champ_Fonc_base)  or OWN_PTR(Champ_Don_base) with P1Bulle discretization" << finl;
      Cerr << "and a multi_scalaire nature." <<finl;
      exit();
    }
}
 
void VEF_discretisation::discretiser_champ_fonc_don_(const Motcle& directive, const Domaine_dis_base& z, Nature_du_champ nature, const Noms& noms, const Noms& unites, int nb_comp, double temps,
                                                     Objet_U& champ) const
{
  // Two pointers for easy access to champ_don or champ_fonc depending on the type of the champ object.
  OWN_PTR(Champ_Fonc_base)  * champ_fonc = dynamic_cast<OWN_PTR(Champ_Fonc_base)*>(&champ);
  OWN_PTR(Champ_Don_base) * champ_don = dynamic_cast<OWN_PTR(Champ_Don_base)*>(&champ);
 
  const Domaine_VEF& domaine_vef = ref_cast(Domaine_VEF, z);
 
  Motcles motcles(8);
  motcles[0] = "pression";    // Standard choice for pressure
  motcles[1] = "temperature"; // Standard choice for temperature
  motcles[2] = "divergence_vitesse"; // Field type obtained by computing div v
  motcles[3] = "champ_elem";  // Create a field at elements (type P0)
  motcles[4] = "vitesse";     // Standard choice for velocity
  motcles[5] = "gradient_pression";  // Field type obtained by computing grad P
  motcles[6] = "champ_sommets"; // Create a field at vertices
  motcles[7] = "champ_face"; // Create a field at faces
 
  // The velocity field type depends on the element type:
  Nom type_champ_vitesse, type_champ_sommets;
  Nom type_elem_domaine = domaine_vef.domaine().type_elem()->que_suis_je();
  if (Motcle(type_elem_domaine) != "Segment")
    {
      const Elem_VEF_base& elem_vef = domaine_vef.type_elem();
      if (sub_type(Tri_VEF, elem_vef) || sub_type(Tetra_VEF, elem_vef))
        {
          type_champ_vitesse = "Champ_Fonc_P1NC";
          type_champ_sommets = "Champ_Fonc_P1_VEF";
        }
      else if (sub_type(Quadri_VEF, elem_vef) || sub_type(Hexa_VEF, elem_vef))
        {
          type_champ_vitesse = "Champ_Fonc_Q1NC";
          type_champ_sommets = "Champ_Fonc_Q1_VEF";
        }
      else
        {
          Cerr << "VEF_discretisation::discretiser_champ :\n The geometric element ";
          Cerr << elem_vef.que_suis_je();
          Cerr << " is not supported." << finl;
          exit();
        }
    }
 
  Nom type;
  int default_nb_comp = 0; // Default number of components
  int rang = motcles.search(directive);
  switch(rang)
    {
    case 0:
      type = "Champ_Fonc_P0_VEF";
      default_nb_comp = 1;
      break;
    case 1:
      type = type_champ_vitesse;
      default_nb_comp = 1;
      break;
    case 2:
      type = "Champ_Fonc_P0_VEF";
      default_nb_comp = 1;
      break;
    case 3:
      type = "Champ_Fonc_P0_VEF";
      default_nb_comp = 1;
      break;
    case 4:
      type = type_champ_vitesse;
      default_nb_comp = dimension;
      break;
    case 5:
      type = type_champ_vitesse;
      default_nb_comp = dimension;
      break;
    case 6:
      type = type_champ_sommets;
      default_nb_comp = 1;
      break;
    case 7:
      type = type_champ_vitesse;
      default_nb_comp = 1;
      break;
 
    default:
      assert(rang < 0);
      break;
    }
 
  if (directive == DEMANDE_DESCRIPTION)
    Cerr << "VEF_discretisation : " << motcles;
 
  // If the directive is not recognized (or if it is a description request)
  // then call the ancestor:
  if (rang < 0)
    {
      if (champ_fonc)
        Discret_Thyd::discretiser_champ(directive, z, nature, noms, unites, nb_comp, temps, *champ_fonc);
      else
        Discret_Thyd::discretiser_champ(directive, z, nature, noms, unites, nb_comp, temps, *champ_don);
      return;
    }
 
  // Compute the number of degrees of freedom
  int nb_ddl = 0;
  if (type == "Champ_Fonc_P0_VEF")
    nb_ddl = z.nb_elem();
  else if (type == type_champ_vitesse)
    nb_ddl = domaine_vef.nb_faces();
  else if ((type == "Champ_Fonc_P1_VEF") || (type == "Champ_Fonc_Q1_VEF"))
    nb_ddl = domaine_vef.nb_som();
  else
    assert(0);
 
  if (nb_comp < 0)
    nb_comp = default_nb_comp;
  assert(nb_comp > 0);
  if (champ_fonc)
    creer_champ(*champ_fonc, z, type, noms[0], unites[0], nb_comp, nb_ddl, temps, directive, que_suis_je());
  else
    creer_champ(*champ_don, z, type, noms[0], unites[0], nb_comp, nb_ddl, temps, directive, que_suis_je());
 
  if ((nature == multi_scalaire) && (champ_fonc))
    {
      champ_fonc->valeur().fixer_nature_du_champ(nature);
      champ_fonc->valeur().fixer_unites(unites);
      champ_fonc->valeur().fixer_noms_compo(noms);
    }
  else if ((nature == multi_scalaire) && (champ_don))
    {
      Cerr << "There is no field of type OWN_PTR(Champ_Don_base) with a multi_scalaire nature." << finl;
      exit();
    }
}
 
void VEF_discretisation::distance_paroi(const Schema_Temps_base& sch, Domaine_dis_base& z, OWN_PTR(Champ_Fonc_base)& ch) const
{
  Cerr << "Discretization of the wall distance" << finl;
  Domaine_VEF& domaine_vef = ref_cast(Domaine_VEF, z);
  ch.typer("Champ_Fonc_P0_VEF");
  Champ_Fonc_P0_VEF& ch_dist_paroi = ref_cast(Champ_Fonc_P0_VEF, ch.valeur());
  ch_dist_paroi.associer_domaine_dis_base(domaine_vef);
  ch_dist_paroi.nommer("distance_paroi");
  ch_dist_paroi.fixer_nb_comp(1);
  ch_dist_paroi.fixer_nb_valeurs_nodales(domaine_vef.nb_elem());
  ch_dist_paroi.fixer_unite("m");
  ch_dist_paroi.changer_temps(sch.temps_courant());
}
 
void VEF_discretisation::distance_paroi_globale(const Schema_Temps_base& sch, Domaine_dis_base& z, OWN_PTR(Champ_Fonc_base)& ch) const
{
  Cerr << "Discretization of global wall distance" << finl;
  Noms noms(1), unites(1);
  noms[0] = Nom("distance_paroi_globale");
  unites[0] = Nom("m");
  discretiser_champ(Motcle("champ_elem"), z, scalaire, noms , unites, 1, 0, ch);
}
 
void VEF_discretisation::vorticite(Domaine_dis_base& z, const Champ_Inc_base& ch_vitesse, OWN_PTR(Champ_Fonc_base)& ch) const
{
  Cerr << "Discretization of vorticity " << finl;
  const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, z);
 
  if (sub_type(Tri_VEF,domaine_VEF.type_elem()) || sub_type(Tetra_VEF, domaine_VEF.type_elem()))
    {
      ch.typer("Rotationnel_Champ_P1NC");
      const Champ_P1NC& vit = ref_cast(Champ_P1NC, ch_vitesse);
      Rotationnel_Champ_P1NC& ch_W = ref_cast(Rotationnel_Champ_P1NC, ch.valeur());
      ch_W.associer_domaine_dis_base(domaine_VEF);
      ch_W.associer_champ(vit);
      ch_W.nommer("vorticite");
      if (dimension == 2)
        ch_W.fixer_nb_comp(1);
      else
        {
          ch_W.fixer_nb_comp(dimension);
          ch_W.fixer_nom_compo(0, "vorticitex");
          ch_W.fixer_nom_compo(1, "vorticitey");
          ch_W.fixer_nom_compo(2, "vorticitez");
        }
      ch_W.fixer_nb_valeurs_nodales(domaine_VEF.nb_elem());
      ch_W.fixer_unite("s-1");
      ch_W.changer_temps(ch_vitesse.temps());
    }
  else if (sub_type(Quadri_VEF,domaine_VEF.type_elem()) || sub_type(Hexa_VEF, domaine_VEF.type_elem()))
    {
      ch.typer("Rotationnel_Champ_Q1NC");
      const Champ_Q1NC& vit = ref_cast(Champ_Q1NC, ch_vitesse);
      Rotationnel_Champ_Q1NC& ch_W = ref_cast(Rotationnel_Champ_Q1NC, ch.valeur());
      ch_W.associer_domaine_dis_base(domaine_VEF);
      ch_W.associer_champ(vit);
      ch_W.nommer("vorticite");
      if (dimension == 2)
        ch_W.fixer_nb_comp(1);
      else
        {
          ch_W.fixer_nb_comp(dimension);
          ch_W.fixer_nom_compo(0, "vorticitex");
          ch_W.fixer_nom_compo(1, "vorticitey");
          ch_W.fixer_nom_compo(2, "vorticitez");
        }
      ch_W.fixer_nb_valeurs_nodales(domaine_VEF.nb_elem());
      ch_W.fixer_unite("s-1");
      ch_W.changer_temps(ch_vitesse.temps());
    }
  else
    {
      Cerr << "Problem in element typing in VEF_discretisation::vorticite" << finl;
      exit();
    }
}
 
void VEF_discretisation::creer_champ_vorticite(const Schema_Temps_base& sch, const Champ_Inc_base& ch_vitesse, OWN_PTR(Champ_Fonc_base)& ch) const
{
  if (sub_type(Champ_P1NC, ch_vitesse))
    {
      ch.typer("Rotationnel_Champ_P1NC");
      const Champ_P1NC& vit = ref_cast(Champ_P1NC, ch_vitesse);
      const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, vit.domaine_dis_base());
      Rotationnel_Champ_P1NC& ch_W = ref_cast(Rotationnel_Champ_P1NC, ch.valeur());
      ch_W.associer_domaine_dis_base(domaine_VEF);
      ch_W.associer_champ(vit);
      ch_W.nommer("vorticite");
      if (dimension == 2)
        ch_W.fixer_nb_comp(1);
      else
        {
          ch_W.fixer_nb_comp(dimension);
          ch_W.fixer_nom_compo(0, "vorticitex");
          ch_W.fixer_nom_compo(1, "vorticitey");
          ch_W.fixer_nom_compo(2, "vorticitez");
        }
      ch_W.fixer_nb_valeurs_nodales(domaine_VEF.nb_elem());
      ch_W.fixer_unite("s-1");
      ch_W.changer_temps(sch.temps_courant());
    }
  else if (sub_type(Champ_Q1NC, ch_vitesse))
    {
      ch.typer("Rotationnel_Champ_Q1NC");
      const Champ_Q1NC& vit = ref_cast(Champ_Q1NC, ch_vitesse);
      const Domaine_VEF& domaine_VEF = ref_cast(Domaine_VEF, vit.domaine_dis_base());
      Rotationnel_Champ_Q1NC& ch_W = ref_cast(Rotationnel_Champ_Q1NC, ch.valeur());
      ch_W.associer_domaine_dis_base(domaine_VEF);
      ch_W.associer_champ(vit);
      ch_W.nommer("vorticite");
      if (dimension == 2)
        ch_W.fixer_nb_comp(1);
      else
        {
          ch_W.fixer_nb_comp(dimension);
          ch_W.fixer_nom_compo(0, "vorticitex");
          ch_W.fixer_nom_compo(1, "vorticitey");
          ch_W.fixer_nom_compo(2, "vorticitez");
        }
      ch_W.fixer_nb_valeurs_nodales(domaine_VEF.nb_elem());
      ch_W.fixer_unite("s-1");
      ch_W.changer_temps(sch.temps_courant());
    }
  else
    {
      Cerr << "Problem in element typing in VEF_discretisation::creer_champ_vorticite" << finl;
      exit();
    }
}
 
/*! @brief Discretizes the incompressible fluid in VEF, specifically K and N.
 *
 * @param (Domaine_dis_base&) domain to discretize
 * @param (Fluide_Ostwald&) fluid to discretize
 * @param (Champ_Inc_base&) velocity
 * @param (Champ_Inc_base&) temperature
 */
void VEF_discretisation::proprietes_physiques_fluide_Ostwald(const Domaine_dis_base& z, Fluide_Ostwald& le_fluide, const Navier_Stokes_std& eqn_hydr, const Champ_Inc_base& ch_temper) const
{
  Cerr << "VEF discretization of the Ostwald fluid" << finl;
  const Domaine_VEF& domaine_vef = ref_cast(Domaine_VEF, z);
  const Champ_Inc_base& ch_vitesse = eqn_hydr.inconnue();
  const Champ_P1NC& vit = ref_cast(Champ_P1NC, ch_vitesse);
 
  Champ_Don_base& mu = le_fluide.viscosite_dynamique();
  //  mu is always a champ_Ostwald_VEF, the following must always be done
  Champ_Ostwald_VEF& ch_mu = ref_cast(Champ_Ostwald_VEF, mu);
  Cerr << "associate domainedisbase VEF" << finl;
  ch_mu.associer_domaine_dis_base(domaine_vef);
  ch_mu.associer_fluide(le_fluide);
  ch_mu.associer_champ(vit);
  ch_mu.associer_eqn(eqn_hydr);
  Cerr << "associations completed: domain dis base, fluid, field VEF" << finl;
  ch_mu.fixer_nb_comp(1);
 
  Cerr << "done fixer_nb_valeurs_nodales" << finl;
  Cerr << "nb_valeurs_nodales VEF = " << domaine_vef.nb_elem() << finl;
  ch_mu.fixer_nb_valeurs_nodales(domaine_vef.nb_elem());
 
  Cerr << "done changer_temps" << finl;
  ch_mu.changer_temps(vit.temps());
  Cerr << "mu VEF is discretized " << finl;
}
 
void VEF_discretisation::critere_Q(const Domaine_dis_base& z, const Domaine_Cl_dis_base& zcl, const Champ_Inc_base& ch_vitesse, OWN_PTR(Champ_Fonc_base)& ch) const
{
  // We pass zcl so there is only one method regardless of the discretization,
  // but it is not used here!!!
  Cerr << "Discretization of Q criterion " << finl;
  const Champ_P1NC& vit = ref_cast(Champ_P1NC, ch_vitesse);
  const Domaine_VEF& domaine_vef = ref_cast(Domaine_VEF, z);
  ch.typer("Critere_Q_Champ_P1NC");
  Critere_Q_Champ_P1NC& ch_cQ = ref_cast(Critere_Q_Champ_P1NC, ch.valeur());
  ch_cQ.associer_domaine_dis_base(domaine_vef);
  ch_cQ.associer_champ(vit);
  ch_cQ.nommer("Critere_Q");
  ch_cQ.fixer_nb_comp(1);
  ch_cQ.fixer_nb_valeurs_nodales(domaine_vef.nb_elem());
  ch_cQ.fixer_unite("s-2");
  ch_cQ.changer_temps(ch_vitesse.temps());
}
 
void VEF_discretisation::y_plus(const Domaine_dis_base& z, const Domaine_Cl_dis_base& zcl, const Champ_Inc_base& ch_vitesse, OWN_PTR(Champ_Fonc_base)& ch) const
{
  Cerr << "Discretization of y_plus" << finl;
  const Champ_P1NC& vit = ref_cast(Champ_P1NC, ch_vitesse);
  const Domaine_VEF& domaine_vef = ref_cast(Domaine_VEF, z);
  const Domaine_Cl_VEF& domaine_cl_vef = ref_cast(Domaine_Cl_VEF, zcl);
  ch.typer("Y_plus_Champ_P1NC");
  Y_plus_Champ_P1NC& ch_yp = ref_cast(Y_plus_Champ_P1NC, ch.valeur());
  ch_yp.associer_domaine_dis_base(domaine_vef);
  ch_yp.associer_domaine_Cl_dis_base(domaine_cl_vef);
  ch_yp.associer_champ(vit);
  ch_yp.nommer("Y_plus");
  ch_yp.fixer_nb_comp(1);
  ch_yp.fixer_nb_valeurs_nodales(domaine_vef.nb_elem());
  ch_yp.fixer_unite("adimensionnel");
  ch_yp.changer_temps(ch_vitesse.temps());
}
 
void VEF_discretisation::t_paroi(const Domaine_dis_base& z,const Domaine_Cl_dis_base& zcl, const Champ_Inc_base& ch_temp, OWN_PTR(Champ_Fonc_base)& ch) const
{
  Cerr << "Discretization of temperature_paroi" << finl;
  const Champ_P1NC& temp = ref_cast(Champ_P1NC, ch_temp);
  const Domaine_VEF& domaine_vef = ref_cast(Domaine_VEF, z);
  const Domaine_Cl_VEF& domaine_cl_vef = ref_cast(Domaine_Cl_VEF, zcl);
  ch.typer("T_paroi_Champ_P1NC");
  T_paroi_Champ_P1NC& ch_tp = ref_cast(T_paroi_Champ_P1NC, ch.valeur());
  ch_tp.associer_domaine_dis_base(domaine_vef);
  ch_tp.associer_domaine_Cl_dis_base(domaine_cl_vef);
  ch_tp.associer_champ(temp);
  ch_tp.nommer("temperature_paroi");
  ch_tp.add_synonymous("wall_temperature");
  ch_tp.fixer_nb_comp(1);
  ch_tp.fixer_nb_valeurs_nodales(domaine_vef.nb_elem());
  ch_tp.fixer_unite("K-C");
  ch_tp.changer_temps(ch_temp.temps());
}
 
void VEF_discretisation::grad_u(const Domaine_dis_base& z, const Domaine_Cl_dis_base& zcl, const Champ_Inc_base& ch_vitesse, OWN_PTR(Champ_Fonc_base)& ch) const
{
 
  const Champ_P1NC& vit = ref_cast(Champ_P1NC, ch_vitesse);
  const Domaine_VEF& domaine_vef = ref_cast(Domaine_VEF, z);
  const int nb_comp = dimension * dimension;
  const Domaine_Cl_VEF& domaine_cl_vef = ref_cast(Domaine_Cl_VEF, zcl);
  ch.typer("grad_u_Champ_P1NC");
  grad_U_Champ_P1NC& ch_grad_u = ref_cast(grad_U_Champ_P1NC, ch.valeur());
  ch_grad_u.fixer_nature_du_champ(vectoriel);
  ch_grad_u.associer_domaine_dis_base(domaine_vef);
  ch_grad_u.associer_domaine_Cl_dis_base(domaine_cl_vef);
  ch_grad_u.associer_champ(vit);
  ch_grad_u.nommer("gradient_vitesse");
  ch_grad_u.fixer_nb_comp(nb_comp);
  ch_grad_u.fixer_nb_valeurs_nodales(domaine_vef.nb_elem());
  ch_grad_u.fixer_unite("s-1");
  ch_grad_u.changer_temps(ch_vitesse.temps());
 
  if (dimension == 2)
    {
      ch_grad_u.fixer_nom_compo(0, "dUdX"); // du/dx
      ch_grad_u.fixer_nom_compo(1, "dUdY"); // du/dy
      ch_grad_u.fixer_nom_compo(2, "dVdX"); // dv/dx
      ch_grad_u.fixer_nom_compo(3, "dVdY"); // dv/dy
    }
  else
    {
      ch_grad_u.fixer_nom_compo(0, "dUdX"); // du/dx
      ch_grad_u.fixer_nom_compo(1, "dUdY"); // du/dy
      ch_grad_u.fixer_nom_compo(2, "dUdZ"); // du/dz
      ch_grad_u.fixer_nom_compo(3, "dVdX"); // dv/dx
      ch_grad_u.fixer_nom_compo(4, "dVdY"); // dv/dy
      ch_grad_u.fixer_nom_compo(5, "dVdZ"); // dv/dz
      ch_grad_u.fixer_nom_compo(6, "dWdX"); // dw/dx
      ch_grad_u.fixer_nom_compo(7, "dWdY"); // dw/dy
      ch_grad_u.fixer_nom_compo(8, "dWdZ"); // dw/dz
    }
 
}
 
void VEF_discretisation::grad_T(const Domaine_dis_base& z, const Domaine_Cl_dis_base& zcl, const Champ_Inc_base& ch_temperature, OWN_PTR(Champ_Fonc_base)& ch) const
{
  Cerr << "Discretization of gradient_temperature" << finl;
 
  const Domaine_VEF& domaine_vef = ref_cast(Domaine_VEF, z);
  creer_champ(ch, domaine_vef, "gradient_temperature_Champ_P1NC", "gradient_temperature", "K/m", dimension, domaine_vef.nb_elem(), ch_temperature.temps());
 
  grad_T_Champ_P1NC& ch_gt = ref_cast(grad_T_Champ_P1NC, ch.valeur());
  const Champ_P1NC& temp = ref_cast(Champ_P1NC, ch_temperature);
  const Domaine_Cl_VEF& domaine_cl_vef = ref_cast(Domaine_Cl_VEF, zcl);
  ch_gt.associer_domaine_Cl_dis_base(domaine_cl_vef);
  ch_gt.associer_champ(temp);
}
 
void VEF_discretisation::h_conv(const Domaine_dis_base& z, const Domaine_Cl_dis_base& zcl, const Champ_Inc_base& ch_temperature, OWN_PTR(Champ_Fonc_base)& ch, Motcle& nom, int temp_ref) const
{
  Cerr << "Discretization of h_conv" << finl;
  const Champ_P1NC& temp = ref_cast(Champ_P1NC, ch_temperature);
  const Domaine_VEF& domaine_vef = ref_cast(Domaine_VEF, z);
  const Domaine_Cl_VEF& domaine_cl_vef = ref_cast(Domaine_Cl_VEF, zcl);
  ch.typer("h_conv_Champ_P1NC");
  h_conv_Champ_P1NC& ch_gt = ref_cast(h_conv_Champ_P1NC, ch.valeur());
  ch_gt.associer_domaine_dis_base(domaine_vef);
  ch_gt.associer_domaine_Cl_dis_base(domaine_cl_vef);
  ch_gt.associer_champ(temp);
  ch_gt.temp_ref() = temp_ref;
  ch_gt.nommer(nom);
  ch_gt.fixer_nb_comp(1);
  ch_gt.fixer_nb_valeurs_nodales(domaine_vef.nb_elem());
  ch_gt.fixer_unite("W/m2.K");
  ch_gt.changer_temps(ch_temperature.temps());
}
 
void VEF_discretisation::taux_cisaillement(const Domaine_dis_base& z, const Domaine_Cl_dis_base& zcl, const Champ_Inc_base& ch_vitesse, OWN_PTR(Champ_Fonc_base)& champ) const
{
  const Domaine_VEF& domaine_vef = ref_cast(Domaine_VEF, z);
  const Domaine_Cl_VEF& domaine_cl_vef = ref_cast(Domaine_Cl_VEF, zcl);
  champ.typer("Taux_cisaillement_P0_VEF");
  Taux_cisaillement_P0_VEF& ch = ref_cast(Taux_cisaillement_P0_VEF, champ.valeur());
  ch.associer_domaine_dis_base(domaine_vef);
  const Champ_P1NC& vit = ref_cast(Champ_P1NC, ch_vitesse);
  ch.associer_champ(vit, domaine_cl_vef);
  ch.nommer("Taux_cisaillement");
  ch.fixer_nb_comp(1);
  ch.fixer_nb_valeurs_nodales(domaine_vef.nb_elem());
  ch.fixer_unite("s-1");
  ch.changer_temps(ch_vitesse.temps());
}
 
void VEF_discretisation::modifier_champ_tabule(const Domaine_dis_base& domaine_dis, Champ_Fonc_Tabule& le_champ_tabule, const VECT(OBS_PTR(Champ_base)) &ch_inc) const
{
  le_champ_tabule.typer_champ_tabule_discretise("Champ_Fonc_Tabule_P0_VEF");
  Champ_Fonc_Tabule_P0_VEF& le_champ_tabule_dis = ref_cast(Champ_Fonc_Tabule_P0_VEF, le_champ_tabule.le_champ_tabule_discretise());
  le_champ_tabule_dis.associer_domaine_dis_base(domaine_dis);
  le_champ_tabule_dis.associer_param(ch_inc, le_champ_tabule.table());
  le_champ_tabule_dis.nommer(le_champ_tabule.le_nom()); // We give a name to this field, help for debug
  le_champ_tabule_dis.fixer_nb_comp(le_champ_tabule.nb_comp());
  le_champ_tabule_dis.fixer_nb_valeurs_nodales(domaine_dis.nb_elem());
  le_champ_tabule_dis.changer_temps(ch_inc[0]->temps());
}
 
void VEF_discretisation::residu(const Domaine_dis_base& z, const Champ_Inc_base& ch_inco, OWN_PTR(Champ_Fonc_base)& champ) const
{
  Nom ch_name(ch_inco.le_nom());
  ch_name += "_residu";
  Cerr << "Discretization of " << ch_name << finl;
 
  const Domaine_VEF& domaine_vef = ref_cast(Domaine_VEF, z);
  int nb_comp = ch_inco.valeurs().line_size();
  Discretisation_base::discretiser_champ("champ_face", domaine_vef, ch_name, "units_not_defined", nb_comp, ch_inco.temps(), champ);
  Champ_Fonc_base& ch_fonc = ref_cast(Champ_Fonc_base, champ.valeur());
  DoubleTab& tab = ch_fonc.valeurs();
  tab = -10000.0;
  Cerr << "[Information] Discretisation_base::residu : the residue is set to -10000.0 at initial time" << finl;
}
 
Domaine_dis_base& VEF_discretisation::discretiser() const
{
  Domaine_dis_base& dom = Discretisation_base::discretiser();
  Domaine_VEF& zvef = ref_cast(Domaine_VEF, dom);
  zvef.discretiser_suite(*this);
  return zvef;
}