Cahn_Hilliard.cpp

/****************************************************************************
* Copyright (c) 2021, 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 <Cahn_Hilliard.h>
#include <Schema_Cahn_Hilliard.h>
#include <Discretisation_base.h>
#include <Domaine_VDF.h>
#include <Discret_Thyd.h>
#include <Probleme_base.h>
#include <Cond_lim_base.h>
#include <Cond_lim.h>
#include <Champ_front_uniforme.h>
#include <CL_Types_include.h>
#include <Courant_impose.h>
#include <EChaine.h>
#include <Avanc.h>
#include <Debog.h>
#include <TRUST_2_PDI.h>
#include <Milieu_Incompressible_Phase_Field.h>
 
Implemente_instanciable(Cahn_Hilliard,"Cahn_Hilliard",Equation_base);
 
Sortie& Cahn_Hilliard::printOn(Sortie& os) const
{
  os << le_dom_Cl_dis_mutilde_ << finl;
  os << mutilde_ << finl;
  if (mobilite_explicite_)
    os << " -- mobilite explicite = TRUE " << finl;
  else
    os << " -- mobilite explicite = FALSE " << finl;
  return Equation_base::printOn(os);
}
 
Entree& Cahn_Hilliard::readOn(Entree& is)
{
  Equation_base::readOn(is);
 
  // Ici on peut rajouter des set_fichier(...) pour sortir les valeurs des termes diffusifs de l'équation
//    terme_mobilite_.set_fichier("Convection_Cahn_Hilliard");
//    terme_mobilite_.set_description("Momentum flow rate [N] if SI units used");
 
  // XXX HYPER IMPORTANT : On fixe set_diffusion_multiscalaire à TRUE
  set_diffusion_multi_scalaire(true);
  Cerr << "EQUATION CAHN HILLIARD : Setting multiscalar boolean to TRUE." << finl;
 
  return is;
 
}
 
void Cahn_Hilliard::set_param(Param& param) const
{
  Equation_base::set_param(param);
  // Spécifiques à l'équation de Cahn-Hilliard :
  param.ajouter_flag("mobilite_explicite",&mobilite_explicite_);
  param.ajouter_flag("passer_par_ajouter_kappa",&passer_par_ajouter_kappa_);
  param.ajouter_flag("passer_par_ajouter_mobilite",&passer_par_ajouter_mobilite_);
  param.ajouter_non_std("conditions_limites_potentiel|boundary_conditions_potential",(this),Param::REQUIRED);
}
 
int Cahn_Hilliard::lire_motcle_non_standard(const Motcle& mot, Entree& is)
{
  if (mot=="conditions_limites_potentiel|boundary_conditions_potential")
    {
      lire_cl_potentiel(is);
      verif_Cl_mutilde();
      return 1;
    }
  else
    return Equation_base::lire_motcle_non_standard(mot,is);
}
 
/*! @brief Lecture des conditions limites sur un flot d'entree.
 *
 * voir Domaine_Cl_dis_base::readOn
 *
 * @param (Entree& is) le flot d'entree
 * @return (Entree&) le flot d'entree modifie
 * @throws le domaine des conditions aux limites discretisee est vide
 */
Entree& Cahn_Hilliard::lire_cl_potentiel(Entree& is)
{
  Cerr << "Reading of boundaries conditions on chemical potential \n";
  if(!le_dom_Cl_dis_mutilde_)
    {
      Cerr << "Error while reading boundaries conditions : " <<
           que_suis_je() << finl;
      Cerr << "The Domaine_Cl_dis_base corresponding to mutilde is empty ..." << finl;
      exit();
    }
  is >> le_dom_Cl_dis_mutilde_.valeur();
  return is;
}
 
 
/*! @brief Renvoie le nombre d'operateurs de l'equation: Pour Cahn-Hilliard c'est 2.
 *
 * @return (int) le nombre d'operateur de l'equation
 */
int Cahn_Hilliard::nombre_d_operateurs() const
{
  return 2;
}
 
int Cahn_Hilliard::nombre_d_operateurs_tot() const
{
  return nombre_d_operateurs();
}
 
/*! @brief Renvoie le i-eme operateur de l'equation:
 *
 * @param (int i) l'index de l'operateur a renvoyer
 * @return (Operateur&) l'operateur indexe par i
 */
const Operateur& Cahn_Hilliard::operateur(int i) const
{
  switch(i)
    {
    case 0:
      return terme_mobilite_;
    case 1:
      return terme_kappa_;
    default :
      Cerr << "Error for Cahn_Hilliard::operateur(int i)" << finl;
      Cerr << "Cahn_Hilliard has " << nombre_d_operateurs() <<" operators "<<finl;
      Cerr << "and you are trying to access the " << i <<" th one."<< finl;
      exit();
    }
  // Pour les compilos!!
  return terme_mobilite_;
}
 
/*! @brief Renvoie le i-eme operateur de l'equation:
 * @param (int i) l'index de l'operateur a renvoyer
 * @return (Operateur&) l'operateur indexe par i
 */
Operateur& Cahn_Hilliard::operateur(int i)
{
  switch(i)
    {
    case 0:
      return terme_mobilite_;
    case 1:
      return terme_kappa_;
    default :
      Cerr << "Error for Cahn_Hilliard::operateur(int i)" << finl;
      Cerr << "Cahn_Hilliard has " << nombre_d_operateurs() <<" operators "<<finl;
      Cerr << "and you are trying to access the " << i <<" th one."<< finl;
      exit();
    }
  // Pour les compilos!!
  return terme_mobilite_;
}
 
/*! @brief Renvoie le nom du domaine d'application: "Cahn_Hilliard".
 *
 * @return (Motcle&) lenom representant le domaine d'application
 */
const Motcle& Cahn_Hilliard::domaine_application() const
{
  static Motcle domaine = "Cahn_Hilliard";
  return domaine;
}
 
/*! @brief Renvoie la concentration (champ inconnue de l'equation) (version const)
 *
 * @return (Champ_Inc_base&) le champ inconnue representant la concentration
 */
const Champ_Inc_base& Cahn_Hilliard::inconnue() const
{
  return concentration_.valeur();
}
 
/*! @brief Renvoie la concentration (champ inconnue de l'equation)
 *
 * @return (Champ_Inc_base&) le champ inconnue representant la concentration
 */
Champ_Inc_base& Cahn_Hilliard::inconnue()
{
  return concentration_.valeur();
}
 
/*! @brief S'associe au probleme: apelle Equation_base::associer_pb_base(const Probleme_base&)
 *  les deux termes diffusifs doivent aussi être initialisés puisqu'ils ne sont pas lus dans le JDD
 * @param (Probleme_base& pb) le probleme auquel s'associer
 */
void Cahn_Hilliard::associer_pb_base(const Probleme_base& pb)
{
  Equation_base::associer_pb_base(pb);
}
 
/*! @brief Associe un mileu physique a l'equation en construisant dynamiquement (cast) un objet de type Milieu_Incompressible_Phase_Field
 *
 *     a partir de l'objet Milieu_base passe en parametre.
 *
 * @param (Milieu_base& un_milieu) le milieu a associer a l'equation
 */
void Cahn_Hilliard::associer_milieu_base(const Milieu_base& un_milieu)
{
  if (sub_type(Milieu_Incompressible_Phase_Field, un_milieu))
    {
      const Milieu_Incompressible_Phase_Field& un_fluide_ch = ref_cast(Milieu_Incompressible_Phase_Field,un_milieu);
      milieu_ = un_fluide_ch;
      fermeture_ = un_fluide_ch.fermeture();
    }
  else
    {
      Cerr << "Error of medium type for the method Cahn_Hilliard::associer_milieu_base." << finl;
      Cerr << "You have to implement a new compatible kind of Milieu_base." << finl;
      exit();
    }
}
 
/*! @brief Renvoie le milieu physique de l'equation
 *
 * @return (Milieu_base&) Le milieu
 */
const Milieu_base& Cahn_Hilliard::milieu() const
{
  if (!milieu_)
    {
      Cerr << "You forgot to associate a phase-field medium to the problem named " << probleme().le_nom() << finl;
      Process::exit();
    }
  return milieu_.valeur();
}
 
/*! @brief Renvoie le milieu physique de l'equation
 *
 * @return (Milieu_base&) Le milieu
 */
Milieu_base& Cahn_Hilliard::milieu()
{
  if (!milieu_)
    {
      Cerr << "You forgot to associate a phase-field medium to the problem named " << probleme().le_nom() << finl;
      Process::exit();
    }
  return milieu_.valeur();
}
 
/*! @brief Dicretise l'equation.
 *
 */
void Cahn_Hilliard::discretiser()
{
  Cerr << "!! WARNING !! : for now, equation discretization (Cahn_Hilliard::discretiser) is only for hydraulic problems." << finl;
 
  // C'est ici qu'est initialisé le nombre de paramètres d'ordre par rapport aux fermeture du système
  nb_parametres_d_ordre_ = fermeture().nb_parametres_d_ordre();
  Cerr << "Number of components for each field = " << nb_parametres_d_ordre_ << finl;
 
  // Discretiser la concentration
  discretiser_concentration();
  champs_compris_.ajoute_champ(concentration_);
 
  // Discretiser mutilde
  discretiser_mutilde();
  champs_compris_.ajoute_champ(mutilde_);
 
  if (!kappa_)
    discretiser_kappa();
  champs_compris_.ajoute_champ(kappa_);
 
  if (!mobilite_)
    discretiser_mobilite();
  champs_compris_.ajoute_champ(mobilite_);
 
  Cerr << "Discretizing of the boundary conditions of mutilde... ";
  le_dom_Cl_dis_mutilde_.typer(discretisation().domaine_cl_dis_type());
  le_dom_Cl_dis_mutilde_->associer(domaine_dis());
  Cerr << " OK " << finl;
 
  le_dom_Cl_dis_mutilde_->associer_eqn(*this);
  le_dom_Cl_dis_mutilde_->associer_inconnue(mutilde_);
 
  Equation_base::discretiser();
 
}
 
void Cahn_Hilliard::discretiser_concentration()
{
  const Discret_Thyd& dis=ref_cast(Discret_Thyd, discretisation());
  // Le nombre de paramètres d'ordre définit le nombre de composantes des champs vectoriels d'inconnues
  Cerr << "Concentration discretization... " << finl;
  // Remarque : je ne passe pas par la méthode "concentration" de Discret_Thyd
  dis.discretiser_champ("champ_elem", domaine_dis(), "concentration", ".", nb_parametres_d_ordre_, schema_temps().nb_valeurs_temporelles(), schema_temps().temps_courant(), concentration_);
  concentration_->nommer("concentration");
  concentration_->fixer_nature_du_champ(multi_scalaire); // Important !
 
}
 
void Cahn_Hilliard::discretiser_mutilde()
{
  const Discret_Thyd& dis=ref_cast(Discret_Thyd, discretisation());
  Cerr << "Mutilde discretization... " << finl;
  // On force mutilde à n'avoir qu'une valeur temporelle (pas besoin de passé et futur)
  // Ça reste un Champ_Inc car il est associé à un opérateur de diffusion
  dis.discretiser_champ("champ_elem", domaine_dis(), "potentiel_chimique_generalise", ".", nb_parametres_d_ordre_, 1, schema_temps().temps_courant(), mutilde_);
  mutilde_->nommer("potentiel_chimique_generalise");
  mutilde_->fixer_nature_du_champ(multi_scalaire); // Important !
 
}
 
void Cahn_Hilliard::discretiser_kappa()
{
  const Discret_Thyd& dis=ref_cast(Discret_Thyd, discretisation());
  Cerr << "Kappa discretization... " << finl;
  // Autre possibilité : "Champ_sommets"
  dis.discretiser_champ("champ_elem",domaine_dis(),"capillarite",".",nb_parametres_d_ordre_*nb_parametres_d_ordre_,0.,kappa_);
  kappa_->nommer("capillarite");
  kappa_->fixer_nature_du_champ(multi_scalaire); // Important !
 
}
 
 
void Cahn_Hilliard::discretiser_mobilite()
{
  const Discret_Thyd& dis=ref_cast(Discret_Thyd, discretisation());
  Cerr << "Mobility discretization... " << finl;
  dis.discretiser_champ("champ_elem",domaine_dis(),"mobilite",".",nb_parametres_d_ordre_*nb_parametres_d_ordre_,0.,mobilite_);
  mobilite_->nommer("mobilite");
  mobilite_->fixer_nature_du_champ(multi_scalaire); // Important !
 
}
 
void Cahn_Hilliard::construire_kappa()
{
  // Cerr << "Kappa field is built..." << finl;
  // Cerr << "Dimension of Kappa = " << kappa_->valeurs().dimension(0) << " x " << kappa_->valeurs().dimension(1) << finl;
  for (int x = 0; x < domaine_dis().nb_elem(); x++)
    for (int j = 0; j < nb_parametres_d_ordre_*nb_parametres_d_ordre_; j++)
      kappa_->valeurs()(x, j) = fermeture().kappa()(j);
 
  kappa_->valeurs().echange_espace_virtuel();
  // Cerr << "Kappa field well computed !" << finl;
}
 
void Cahn_Hilliard::construire_mobilite(const DoubleTab& c)
{
  DoubleTrav result_mob = fermeture().mobilite(c);
  result_mob = fermeture().mobilite(c);
  for (int x = 0; x < domaine_dis().nb_elem(); x++)
    for (int j = 0; j < nb_parametres_d_ordre_*nb_parametres_d_ordre_; j++)
      mobilite_->valeurs()(x,j) = result_mob(x,j);
 
  mobilite_->valeurs().echange_espace_virtuel();
}
 
void Cahn_Hilliard::lire_terme_mobilite()
{
  // Opérateur 1 : mobilité
  terme_mobilite_.associer_eqn(*this);
  terme_mobilite_.associer_diffusivite(mobilite_);
 
  Cerr << "[COMPLETER()] Typing of the diffusion operator of mobility... " << finl;
  // Manière adhoc de lire le type d'opérateur de diffusion (soit { } soit { negligeable }
  EChaine diff("{ }");
  diff >> terme_mobilite_;
 
  terme_mobilite_.associer_diffusivite_pour_pas_de_temps(mobilite_);
  // XXX Luis : ATTENTION, big probleme quand on associe un champ avec un nom différent de l'inconnue de l'équation
  // Plus tard il faut utiliser ajouter_blocs et non pas contribuer_a_avec, ni calculer
  terme_mobilite_.associer_champ(mutilde_,"potentiel_chimique_generalise");
 
}
 
void Cahn_Hilliard::lire_terme_kappa()
{
  // Opérateur 2 : capillarité
  terme_kappa_.associer_eqn(*this);
  terme_kappa_.associer_diffusivite(kappa_);
 
  Cerr << "[COMPLETER()] Typing of the diffusion operator of capillarity... " << finl;
  EChaine diff2("{ }");
  diff2 >> terme_kappa_;
 
  terme_kappa_.associer_diffusivite_pour_pas_de_temps(kappa_);
  terme_kappa_.associer_champ(concentration_,"concentration");
}
 
/*! @brief Complete l'equation base, associe le potentiel chimique généralisé a l'equation,
 *     initialise le vecteur résidu et construit les opérateurs de diffusion.
 */
void Cahn_Hilliard::completer()
{
  // On bypass la méthode completer de Equation_base
  Cerr << "####### CAHN_HILLIARD::COMPLETER() ########" << finl;
 
  inconnue().associer_eqn(*this);
  inconnue().completer(le_dom_Cl_dis.valeur());
 
  // CL sur la concentration
  if (le_dom_Cl_dis)
    le_dom_Cl_dis->completer();
 
  mutilde_->associer_eqn(*this);
  mutilde_->completer(le_dom_Cl_dis_mutilde_.valeur());
 
  // CL sur mutilde
  if (le_dom_Cl_dis_mutilde_)
    le_dom_Cl_dis_mutilde_->completer();
 
  // Association du champ d'inconnue (concentration) à son domaine
  inconnue().associer_domaine_cl_dis(le_dom_Cl_dis);
  Cerr << "[COMPLETER()] The unknown field c is associated to boundary conditions." << finl;
 
  // Association du champ mutilde à son domaine
  mutilde_->associer_domaine_cl_dis(le_dom_Cl_dis_mutilde_);
  Cerr << "[COMPLETER()] The potential field mutilde is associated to boundary conditions." << finl;
 
  if (mon_probleme->is_coupled())
    {
      Cerr << "Cahn_Hilliard :: Coupled problem not coded..." << finl;
      exit();
    }
 
  //Ajout d un element vide qui sera renvoye si pas de modele trouve
  RefObjU mod;
  if (liste_modeles_.size()==0)
    liste_modeles_.add(mod);
 
 
  Debog::verifier("Cahn_Hilliard::completer() avant concentration ", concentration_->valeurs());
  Debog::verifier("Cahn_Hilliard::completer() avant mobilite ", mobilite_->valeurs());
  Debog::verifier("Cahn_Hilliard::completer() avant kappa ", kappa_->valeurs());
 
  Cerr << "Mobility field is built..." << finl;
  construire_mobilite(concentration_->valeurs()); // Construit le CHAMP de mobilité
  Cerr << "Mobility field well computed !" << finl;
  Cerr << "Kappa field is built..." << finl;
  construire_kappa(); // Construit le CHAMP de kappa
  Cerr << "Kappa field well computed !" << finl;
 
  lire_terme_mobilite();
  lire_terme_kappa();
 
  Debog::verifier("Cahn_Hilliard::completer() apres concentration ", concentration_->valeurs());
  Debog::verifier("Cahn_Hilliard::completer() apres mobilite ", mobilite_->valeurs());
  Debog::verifier("Cahn_Hilliard::completer() apres kappa ", kappa_->valeurs());
 
 
  // CE QU'IL Y A DANS COMPLETER DE OPERATEUR_BASE
  // On decline le residu par composante, qui correspondent au nombre de paramètres d'ordre (et donc d'équations de Cahn-Hilliard)
  initialise_residu();
  int nb_op=nombre_d_operateurs();
  Cerr << "Number of operators = " << nb_op << finl;
  Nom msg="Equation_base::completer(), nb_op = " ;
  Debog::verifier(msg, nb_op);
  for(int i=0; i<nb_op; i++)
    operateur(i).completer();
 
  if (solveur_masse)
    solv_masse().completer();
 
  les_sources.completer();
  schema_temps().completer();
 
//  Schema_Cahn_Hilliard schee = ref_cast(Schema_Cahn_Hilliard,schema_temps());
//  if (abs(schee.coeff_theta()) < 1.e-12)
  implicite_=0;
//  else
//    implicite_=1;
 
}
 
/*! @brief cf Equation_base::preparer_calcul()
 *
 * @return (int) renvoie toujours 1
 */
int Cahn_Hilliard::preparer_calcul()
{
  const double temps = schema_temps().temps_courant();
  Equation_base::preparer_calcul();
  sources().mettre_a_jour(temps);
 
  mutilde().verifie_valeurs_cl();
  mutilde().changer_temps(temps);
  le_dom_Cl_dis_mutilde_->initialiser(temps);
 
  domaine_Cl_dis_mutilde().imposer_cond_lim(mutilde(),temps);
  mutilde().valeurs().echange_espace_virtuel();
 
  le_dom_Cl_dis->mettre_a_jour(temps);
  le_dom_Cl_dis_mutilde_->mettre_a_jour(temps);
 
  // Mise a jour mutilde
  mutilde_->changer_temps(temps);
 
  Debog::verifier("Cahn_Hilliard::preparer_calcul, concentration", inconnue());
  Debog::verifier("Cahn_Hilliard::preparer_calcul, mutilde", mutilde());
 
  return 1;
}
 
void Cahn_Hilliard::abortTimeStep()
{
  mutilde().abortTimeStep();
  Equation_base::abortTimeStep();
}
 
/* @brief Override. Reset mutilde too !
 */
void Cahn_Hilliard::resetTime(double time)
{
  mutilde().resetTime(time);
  domaine_Cl_dis_mutilde().resetTime(time);
  Equation_base::resetTime(time);
}
 
bool Cahn_Hilliard::initTimeStep(double dt)
{
 
  const Schema_Temps_base& sch_tps = le_schema_en_temps.valeur();
 
  for (int i=1; i<=sch_tps.nb_valeurs_futures(); i++)
    if (i <= mutilde().nb_valeurs_temporelles())
      {
        double tps=sch_tps.temps_futur(i);
        // Mise a jour du temps dans les champs de potentiel chimique
        mutilde().changer_temps_futur(tps,i);
        mutilde().futur(i)=mutilde().valeurs();
        // Mise a jour du temps dans les CL
        domaine_Cl_dis_mutilde().changer_temps_futur(tps,i);
      }
 
  // Mise a jour du temps par defaut des CLs
  domaine_Cl_dis_mutilde().set_temps_defaut(sch_tps.temps_defaut());
 
  bool ddt = Equation_base::initTimeStep(dt);
 
  return ddt;
}
 
/*! @brief cf Equation_base::updateGivenFields()
 *
 * On rajoute ici les mises à jour des CL sur mutilde. Utile ? Je ne sais pas.
 */
bool Cahn_Hilliard::updateGivenFields()
{
  Schema_Temps_base& sch=schema_temps();
  double temps_present=sch.temps_courant();
  double temps_futur=temps_present+sch.pas_de_temps();
 
  // Pour chaque temps futur
  for (int i=1; i<=sch.nb_valeurs_futures(); i++)
    {
      double tps=sch.temps_futur(i);
      // Calcul des CLs a ce temps.
      domaine_Cl_dis_mutilde().mettre_a_jour(tps);
    }
  // Calcul du taux d'accroissement des CLs entre les temps present et futur.
  domaine_Cl_dis_mutilde().calculer_derivee_en_temps(temps_present,temps_futur);
 
  return Equation_base::updateGivenFields();
}
 
/*! @brief for PDI IO: retrieve name, type and dimensions of the fields to save/restore
 *
 */
std::vector<YAML_data> Cahn_Hilliard::data_a_sauvegarder() const
{
  std::vector<YAML_data> data = Equation_base::data_a_sauvegarder();
  std::vector<YAML_data> mutilde_sauv = mutilde_->data_a_sauvegarder();
  data.insert(data.end(), mutilde_sauv.begin(), mutilde_sauv.end());
  return data;
}
 
/*! @brief Appelle Equation_base::sauvegarder(Sortie&) et sauvegarde le potentiel chimique sur un flot de sortie.
 *
 * @param (Sortie& os) un flot de sortie sur lequel sauvegarder
 * @return (int) renvoie toujours 1
 */
int Cahn_Hilliard::sauvegarder(Sortie& os) const
{
  int bytes=0;
  bytes += Equation_base::sauvegarder(os);
  bytes += mutilde_->sauvegarder(os);
  return bytes;
}
 
/*! @brief Effectue une reprise a partir d'un flot d'entree.
 *
 * Appelle Equation_base::reprendre()
 *      et reprend mutilde.
 *
 * @param (Entree& is) un flot d'entree
 * @return (int) renvoie toujours 1
 * @throws la reprise a echoue, identificateur du potentiel chimique non trouve
 */
int Cahn_Hilliard::reprendre(Entree& is)
{
  Equation_base::reprendre(is);
  if(!TRUST_2_PDI::is_PDI_restart())
    {
      double temps = schema_temps().temps_courant();
      Nom ident_mutilde(mutilde_->le_nom());
      ident_mutilde += mutilde_->que_suis_je();
      ident_mutilde += probleme().domaine().le_nom();
      ident_mutilde += Nom(temps,probleme().reprise_format_temps());
      avancer_fichier(is, ident_mutilde);
    }
  mutilde_->reprendre(is);
 
  return 1;
}
 
/*! @brief S'associe au schema_en_temps : peut seulement être Schema_CH_Semi_Implicite
 *  pour l'instant.
 * @param (Schema_Temps_base& un_schema_en_temps) le schema en temps a associer a l'equation
 */
void Cahn_Hilliard::associer_sch_tps_base(const Schema_Temps_base& un_schema_en_temps)
{
  if (sub_type(Schema_Cahn_Hilliard,un_schema_en_temps))
    {
      le_schema_en_temps=un_schema_en_temps;
    }
  else
    {
      Cerr << "Error : " << que_suis_je()
           << " cannot support other time scheme than Schema_Cahn_Hilliard." << finl;
      exit();
    }
}
 
/*! @brief Effectue une mise a jour en temps de l'equation.
 *
 * Appelle Equation_base::mettre_a_jour(double)
 *      et met a jour le potentiel chimique généralisé.
 *
 * @param (double temps) le temps de mise a jour
 */
void Cahn_Hilliard::mettre_a_jour(double temps)
{
  // Mise a jour de la classe mere (on tourne la roue).
  Equation_base::mettre_a_jour(temps);
 
  // Mise a jour de mutilde
  mutilde_->mettre_a_jour(temps);
  // Mise a jour conditions limites de mutilde
  if (temps > schema_temps().temps_courant()) domaine_Cl_dis_mutilde().avancer(temps);
}
 
/*! @brief Calcul du prochain pas de temps.
 *
 *  Luis : Pour l'instant, dt = dt_min, je ne sais pas quel est le pas de temps optimal pour l'équation de Cahn-Hilliard
 *
 * @return (double)
 */
double Cahn_Hilliard::calculer_pas_de_temps() const
{
  // On peut chercher à augmenter le pas de temps en fonction de la norme |c^n+1 - c^n| qui indique si on atteint un niveau stationnaire
  // Mais attention, il peut se passer des mécanismes rapides même après une évolution très lente...
  // On reste sur dt_max pour l'instant
  return le_schema_en_temps->pas_temps_max();
}
 
/*! @brief Returns the time derivative of the unknown I of the equation: dI/dt = M-1*(sum(operators(I) + sources))
 *
 *  In : derivee contains I (but immediatly set to 0)
 *  Out: derivee contains dI/dt
 *
 */
DoubleTab& Cahn_Hilliard::derivee_en_temps_inco(DoubleTab& derivee)
{
  derivee = 0.;
 
  DoubleTrav secmem(derivee);
  // secmem = sum(operators) + sources + equation specific terms
 
  if (implicite_==0)
    {
      // On cherche à construire : dc/dt = M-1 * (D mutilde^n)
      // avec mutilde^n = beta*dW/dc^n - Kc^n
      // où D est la matrice de l'opérateur lié à mu (terme_mobilite_)
      // et K est la matirce de l'opérateur lié à kappa (terme_kappa_)
 
      // On calcule d'abord mutilde^n = beta*dW/dc^n - Kc^n (fonction compute_mutilde)
      compute_mutilde();
      // Ici normalement les valeurs de mutilde (qui est l'inconnue de l'opérateur 0) sont mises à jour
      operateur(0).calculer(secmem); // secmem = D mutilde^n
      les_sources.ajouter(secmem); // secmem += sources
      solveur_masse->appliquer(secmem); // M-1 * secmem
      derivee = secmem;
      derivee.echange_espace_virtuel();
 
    }
  else if (implicite_>0)
    {
//      // En implicite
//      // M dc/dt + D*mutilde^# = f;
//      // où # désigne le temps futur (souvent n+1)
//      // On va d'abord construire la fonction de forme :
//      // H(c) = c^# - c^n - dt * M_c^-1 * D * (M_mu^-1 K c^# + dW/dc^#)
//      // Puis résoudre le problème non linéaire par méthode de Newton : H(c^#) = 0
//      // Enfin, dc/dt = (c^# - c^n) / dt
//      double dt=schema_temps().pas_de_temps();
//
//      // On stocke les valeurs du temps présent dans le tableau du temps passé
//      concentration_->passe() = concentration_->valeurs();
//      compute_mutilde();
//      mutilde_->passe() = mutilde_->valeurs();
//
//      derivee-=concentration_->passe();
//      derivee/=dt;
      Cerr << "[Cahn_Hilliard] Warning: Be careful with implicit part in derivative computation." << finl;
      Process::exit();
    }
  else
    {
      Cerr << "Error in Cahn_Hilliard::derivee_en_temps_inco" << finl;
      Cerr << "implicite_ = " << implicite_ << " has not been initialized!" << finl;
      Cerr << "May be " << que_suis_je() << "::completer() method doesn't call Cahn_Hilliard::completer()" << finl;
      Process::exit();
    }
 
  return derivee;
}
 
/*! @brief Calcul de mutilde = beta*dW/dc - ∇. kappa ∇c
   */
void Cahn_Hilliard::compute_mutilde()
{
  Debog::verifier("Cahn_Hilliard::compute_mutilde() avant mutilde_ ", mutilde_->valeurs());
  Debog::verifier("Cahn_Hilliard::compute_mutilde() avant concentration_ ", concentration_->valeurs());
 
  DoubleTrav terme_non_lin(concentration_->valeurs());
 
  const int nb_elem = terme_non_lin.dimension(0); // reel
  const int nb_param = terme_non_lin.line_size();
 
  if (passer_par_ajouter_kappa_)
    {
      // operateur(1) = l'opérateur de diffusion lié à kappa
      // calculer renvoie : mutilde = CL(c) - K*c
      operateur(1).calculer(mutilde_->valeurs());
 
      // Division par le volume des voxels
      tab_divide_any_shape(mutilde_->valeurs(), ref_cast(Domaine_VF,domaine_dis()).volumes());
 
      // On ajoute la partie liée au potentiel chimique
      terme_non_lin = fermeture().potentiel_chimique().dWdc(concentration_->valeurs());
 
      for (int elem = 0; elem < nb_elem; elem++)
        for (int p = 0; p < nb_param; p++)
          {
            terme_non_lin(elem,p) *= fermeture().beta()(p);
            mutilde_->valeurs()(elem,p) = terme_non_lin(elem,p) - mutilde_->valeurs()(elem,p);
          }
    }
  else // On ne passe pas par la méthode ajouter() de l'opérateur, mais par un produit matrice-vecteur
    {
      DoubleTrav second_terme(concentration_->valeurs());
 
      terme_non_lin = fermeture().potentiel_chimique().dWdc(concentration_->valeurs());
 
      // Calcul du produit matrice / vecteur utilise
      mat_kappa_.multvect_(concentration_->valeurs(), second_terme);
 
      for (int elem = 0; elem < nb_elem; elem++)
        for (int p = 0; p < nb_param; p++)
          {
            terme_non_lin(elem,p) *= fermeture().beta()(p);
            // Attention au signe devant le terme K.c
            mutilde_->valeurs()(elem,p) = terme_non_lin(elem,p) + second_terme(elem,p) - contribution_CL_concentration_(elem,p);
          }
    }
 
  mutilde_->valeurs().echange_espace_virtuel();
  Debog::verifier("Cahn_Hilliard::compute_mutilde() apres mutilde_ ", mutilde_->valeurs());
  Debog::verifier("Cahn_Hilliard::compute_mutilde() apres concentration_ ", concentration_->valeurs());
}
 
 
/*! @brief Calcul de mutilde = beta*dW/dc - ∇. kappa ∇c
 * Stockage dans un autre vecteur
   */
DoubleTab Cahn_Hilliard::compute_mutilde(const DoubleTab& c)
{
  Debog::verifier("Cahn_Hilliard::compute_mutilde_tab() avant mutilde_ ", mutilde_->valeurs());
  DoubleTrav terme_non_lin(concentration_->valeurs());
 
  const int nb_elem = terme_non_lin.dimension(0); // reel
  const int nb_param = terme_non_lin.line_size();
 
  DoubleTrav muti(concentration_->valeurs());
 
  if (passer_par_ajouter_kappa_)
    {
      // operateur(1) = l'opérateur de diffusion lié à kappa
      // calculer renvoie : mutilde = CL(c) - K*c
      operateur(1).calculer(c,muti);
 
      // Division par le volume des voxels
      tab_divide_any_shape(muti, ref_cast(Domaine_VF,domaine_dis()).volumes());
 
      // On ajoute la partie liée au potentiel chimique
      terme_non_lin = fermeture().potentiel_chimique().dWdc(c);
 
      for (int elem = 0; elem < nb_elem; elem++)
        for (int p = 0; p < nb_param; p++)
          {
            terme_non_lin(elem,p) *= fermeture().beta()(p);
            muti(elem,p) = terme_non_lin(elem,p) - muti(elem,p);
          }
    }
  else // On ne passe pas par la méthode ajouter() de l'opérateur, mais par un produit matrice-vecteur
    {
      DoubleTrav second_terme(concentration_->valeurs()) ;
 
      terme_non_lin = fermeture().potentiel_chimique().dWdc(c);
 
      // Calcul du produit matrice / vecteur utilise
      mat_kappa_.multvect_(c, second_terme);
 
      for (int elem = 0; elem < nb_elem; elem++)
        for (int p = 0; p < nb_param; p++)
          {
            terme_non_lin(elem,p) *= fermeture().beta()(p);
            // Attention au signe devant le terme K.c
            muti(elem,p) = terme_non_lin(elem,p) + second_terme(elem,p) - contribution_CL_concentration_(elem,p);
          }
    }
 
  muti.echange_espace_virtuel();
  Debog::verifier("Cahn_Hilliard::compute_mutilde_tab() apres mutilde_ ", mutilde_->valeurs());
  return muti;
}
 
/*! @brief Construit la fonction résidu pour un algorithme de Newton :
 *
 *  En schéma Euler implicite
 *  c^n+1 - c^n - Δt M^-1 D mu^n+1 + CL(mu)^n+1 = 0
 *  avec mu^n+1 = beta*dW/dc^n+1 - M^-1 K c^n+1 - CL(c)^n+1
 *
 *  Nécessite d'être généralisé à toute équation et tout schéma en temps
 *
 */
DoubleTab Cahn_Hilliard::fonction_residu(const DoubleTab& c)
{
  // Variables temporelles
  double dt=schema_temps().pas_de_temps();
  double theta=ref_cast(Schema_Cahn_Hilliard,schema_temps()).coeff_theta();
 
  const int nb_elem = domaine_dis().nb_elem(); // reel
  const int nb_param = nb_parametres_d_ordre();
 
  const DoubleTab& c_n = inconnue().passe(); // Concentration au temps n
 
  DoubleTrav mutilde_theta(concentration_->valeurs());
  DoubleTrav residu(concentration_->valeurs());
 
  if (!mobilite_explicite_)
    update_terme_mobilite(c);
 
  // On calcule d'abord mutilde^theta = beta*dW/dc^n - Kc^n
  // Normalement, mutilde est calculé par rapport aux valeurs de inconnues().valeurs()
  mutilde_theta = compute_mutilde(c);
//  Cerr << "mutilde = " << mutilde_theta << finl;
 
  Debog::verifier("Cahn_Hilliard::fonction_residu() avant c ", c);
 
  // TRES IMPORTANT : La contribution des conditions limites sur mutilde doit être
  // calculée par rapport au mutilde de l'itération de Newton PRÉCÉDENTE !!!
  // C'est dans le cas où les CL de mutilde dépendent de mutilde ou de c (Courant_impose)
 
  // Vraiment ? Je ne suis pas sûr à 100% de ce que j'écris ici.
  // mu(c) est une fonction de c. On le met à jour avec le c de l'itération de Newton précédente.
  for (int i=0; i<le_dom_Cl_dis_mutilde_->nb_cond_lim(); i++)
    {
      Cond_lim& la_cl=le_dom_Cl_dis_mutilde_->les_conditions_limites(i);
      if (sub_type(Courant_impose,la_cl.valeur()))
        {
//          Cerr << "On rentre dans la condition sur la CL et on met à jour le courant !" << finl;
          ref_cast(Courant_impose, la_cl.valeur()).mettre_a_jour_courant(c, mutilde_theta);
          calculer_contribution_CL_mutilde(c);
        }
    }
 
  if (passer_par_ajouter_mobilite_)
    {
      // Ici normalement l'opérateur op(0) lié à la mobilité doit s'appliquer à mutilde_theta
      // operateur(0).calculer(mutilde_theta,tab_residu); // residu = D mutilde^theta
      // Ce n'est pas ce qui est fait !!
      // Car le nom de l'inconnue de l'équation n'est pas le même que celui de l'opérateur.
      // On bypass la méthode surchargée calculer() et on appelle directement ajouter_bloc avec le bon tag.
      terme_mobilite_->ajouter_blocs( { }, residu, { { "potentiel_chimique_generalise", mutilde_theta } });
 
      solv_masse().appliquer(residu); // M-1 * residu
 
      for (int elem = 0; elem < nb_elem; elem++)
        for (int p = 0; p < nb_param; p++)
          residu(elem, p) = c(elem, p) - c_n(elem, p) - theta * dt * residu(elem, p);
    }
  else
    {
      matrice_mobilite().multvect_(mutilde_theta, residu);
 
      for (int elem = 0; elem < nb_elem; elem++)
        for (int p = 0; p < nb_param; p++)
          {
            residu(elem, p) = c(elem, p) - c_n(elem, p) + theta * dt * (residu(elem, p) - contribution_CL_mutilde_(elem, p));
            // Attention au signe devant le terme D.mu (:= intermediaire)
          }
    }
 
  residu.echange_espace_virtuel();
 
  Debog::verifier("Cahn_Hilliard::fonction_residu() apres c ", c);
  Debog::verifier("Cahn_Hilliard::fonction_residu() apres residu ", residu);
 
  return residu;
}
 
/*! @brief Renvoie la contribution des conditions limites sur la concentration
 *              sous forme d'un vecteur (multi-composantes)
 */
void Cahn_Hilliard::calculer_contribution_CL_concentration(const DoubleTab& inco)
{
  DoubleTrav temp(inco), temp2(inco); // inco ne sert qu'à dimensionner
 
  if (!matrices_initialisees_)
    {
      Cerr << "[Cahn_Hilliard] You cannot add a BC contribution without calling the initialisation of matrices before." << finl;
      exit();
    }
 
  // C'est tordu de récupérer les CL comme ça... Beaucoup d'opérations
  terme_kappa_.calculer(temp2,temp); // Calcule b = [CL]^n - K.c^n
  tab_divide_any_shape(temp, ref_cast(Domaine_VF,domaine_dis()).volumes()); // b = M^-1 ([CL]^n - K.c^n)
 
  contribution_CL_concentration_ = temp;
  contribution_CL_concentration_.echange_espace_virtuel();
}
 
/*! @brief Renvoie la contribution des conditions limites sur le potentiel chimique
 *              sous forme d'un vecteur (multi-composantes)
 */
void Cahn_Hilliard::calculer_contribution_CL_mutilde(const DoubleTab& inco)
{
  DoubleTrav temp(inco), temp2(inco); // inco ne sert qu'à dimensionner
 
  if (!matrices_initialisees_)
    {
      Cerr << "[Cahn_Hilliard] You cannot add a BC contribution without calling the initialisation of matrices before." << finl;
      exit();
    }
 
  // C'est tordu de récupérer les CL comme ça... Beaucoup d'opérations
  terme_mobilite_.valeur().ajouter_blocs({}, temp, {{ "potentiel_chimique_generalise",temp2 }} );
  solveur_masse->appliquer(temp); // b = M^-1 ([CL]^n - D.mu^n)
 
  contribution_CL_mutilde_ = temp;
  contribution_CL_mutilde_.echange_espace_virtuel();
}
 
void Cahn_Hilliard::initialiser_matrices()
{
  const Domaine_VF& zvf = ref_cast(Domaine_VF,domaine_dis());
  const DoubleVect& vol = zvf.volumes();
  DoubleTrav inv_vol(vol.size_totale()*nb_parametres_d_ordre_);
 
//  Cerr << "Domain discretization: dx = " << ref_cast(Domaine_VDF,zvf).h_x() << " ; dy = " << ref_cast(Domaine_VDF,zvf).h_y() << finl;
//  Cerr << "Voxel volume V = " << vol(0) << finl;
 
  // Définir un tableau temporaire pour construire la matrice (équivalent de contribuer_a_avec)
  DoubleTrav temp(concentration_->valeurs());
  temp = concentration_->valeurs();
 
  Cerr << " . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ." << finl;
  Cerr << ". . . . . . . . . . . INITIALIZATION: Mobility matrix & Kappa matrix . . . . . . . . . ." << finl;
  Cerr << " . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ." << finl;
 
  terme_kappa_.associer_diffusivite(kappa_);
  terme_kappa_->associer_diffusivite(kappa_);
  terme_kappa_->dimensionner_blocs( { { "concentration", &mat_kappa_ } });
  Cerr << "Dimensions of matrix K = " << mat_kappa_.nb_lignes() << " x " << mat_kappa_.nb_colonnes() << finl;
 
  mat_kappa_.sort_stencil(); // XXX sort tab2 ... au cas ou ...
  terme_kappa_->ajouter_blocs( { { "concentration", &mat_kappa_ } }, temp, { });
 
  // Initialisation des valeurs de mutilde à 0
  mutilde_->valeurs() = 0.;
  temp = mutilde_->valeurs();
 
  terme_mobilite_.associer_diffusivite(mobilite_);
  terme_mobilite_->associer_diffusivite(mobilite_);
  terme_mobilite_->dimensionner_blocs( { { "concentration", &mat_mobilite_ } }); // Attention ! On dimensionne pareil que la concentration
  // mais la matrice mat_mobilite_ doit être associée à mutilde
  Cerr << "Dimensions of matrix M  = " << mat_mobilite_.nb_lignes() << " x " << mat_mobilite_.nb_colonnes() << finl;
 
  mat_mobilite_.sort_stencil(); // XXX sort tab2 ... au cas ou ...
  terme_mobilite_->ajouter_blocs( { { "potentiel_chimique_generalise", &mat_mobilite_ } }, temp, { });
 
  // Stockage de l'inverse des volumes de chaque élément de domaine
  for (int i=0; i<vol.size_totale(); i++)
    for (int p=0; p<nb_parametres_d_ordre_; p++)
      inv_vol(p+i*nb_parametres_d_ordre_) = 1./vol(i);
 
  // Division par les volumes
  mat_mobilite_.diagmulmat(inv_vol);
 
//  Cerr << "Matrice M = ";
//  mat_mobilite_.imprimer_formatte(Cerr);
 
  // Division par les volumes
  mat_kappa_.diagmulmat(inv_vol);
 
//  Cerr << "Matrice K = ";
//  mat_kappa_.imprimer_formatte(Cerr);
 
  matrices_initialisees_=true;
  // Initialisation taille des tableaux des contributions
  contribution_CL_concentration_ = concentration_->valeurs();
  contribution_CL_mutilde_ = mutilde_->valeurs();
 
  // Initialisation des contributions aux CLs
  calculer_contribution_CL_concentration(concentration_->valeurs());
  calculer_contribution_CL_mutilde(mutilde_->valeurs());
}
 
void Cahn_Hilliard::update_terme_mobilite(const DoubleTab& c)
{
  construire_mobilite(c); // Construit le champ de mobilité en fonction d'un champ de concentration c (peut être c^n ou c^theta)
  terme_mobilite_.associer_diffusivite(mobilite_);
  terme_mobilite_->associer_diffusivite(mobilite_);
 
  if (!passer_par_ajouter_mobilite_)
    {
      const Domaine_VF& zvf = ref_cast(Domaine_VF,domaine_dis());
      const DoubleVect& vol = zvf.volumes();
      DoubleTrav inv_vol(vol.size_totale()*nb_parametres_d_ordre_);
 
      // IMPORTANT : Remettre la matrice à 0
      mat_mobilite_.clean();
 
      DoubleTrav temp(mutilde_->valeurs());
      // Met à jours la matrice mat_mobilite_
      mat_mobilite_.sort_stencil(); // XXX sort tab2 ... au cas ou ...
      terme_mobilite_->ajouter_blocs({ { "potentiel_chimique_generalise", &mat_mobilite_ } }, temp, { });
 
      // Stockage de l'inverse des volumes de chaque élément de domaine
      for (int i=0; i<vol.size_totale(); i++)
        for (int p=0; p<nb_parametres_d_ordre_; p++)
          inv_vol(p+i*nb_parametres_d_ordre_) = 1./vol(i);
 
      // Division par les volumes
      mat_mobilite_.diagmulmat(inv_vol);
    }
}
 
/*! @brief Méthode pour vérifier la compatibilité entre les CL de la concentration
 * et les CL du potentiel
 */
int Cahn_Hilliard::verif_Cl_mutilde() const
{
  // D'abord check si les CL de mutilde sont compatibles avec l'équation de Cahn-Hilliard
  le_dom_Cl_dis_mutilde_->les_conditions_limites().compatible_avec_eqn(*this);
 
  // A ce stade, les CL de la concentration sont connues, et celles de mutilde aussi.
  // Test de compatibilité inter-conditions (entre c et mutilde).
  // Les seules CL acceptées pour cette équations pour l'instant sont :
  // 1) Periodique
  // 2) Scalaire_impose_paroi (Dirichlet)
  // 3) Neumann_paroi (flux imposé)
 
  Cerr << "** Checking compatibility between boundary conditions on concentration and potential... **" << finl;
 
  // Boucle sur les frontières du domaine
  for (int n_cl = 0 ; n_cl < domaine_dis().domaine().nb_front_Cl() ; n_cl++)
    {
      const Cond_lim& cl_c=le_dom_Cl_dis->les_conditions_limites(n_cl);
      Nom nom_frontiere = domaine_dis().domaine().bord(n_cl).le_nom();
      const Cond_lim& cl_mutilde=le_dom_Cl_dis_mutilde_->les_conditions_limites(n_cl);
 
      // 1) c periodique, mutilde periodique
      if (sub_type(Periodique,cl_c.valeur()) && sub_type(Periodique,cl_mutilde.valeur()))
        {
          Cerr << "---------↧      " << finl;
          Cerr << "     [" << nom_frontiere << "] Boundary conditions on c and mutilde are compatible: both are PERIODIC." << finl;
          Cerr << "         ↪" << finl;
        }
      // 2) c Neumann
      else if (sub_type(Neumann_paroi,cl_c.valeur()))
        {
          Cerr << "---------↧      " << finl;
          if (sub_type(Neumann_paroi,cl_mutilde.valeur())) // 2.1) mutilde Neumann
            Cerr << "     [" << nom_frontiere << "] Boundary conditions on c and mutilde are compatible: both are NEUMANN." << finl;
          else if (sub_type(Scalaire_impose_paroi,cl_mutilde.valeur())) // 2.2) mutilde Dirichlet
            Cerr << "     [" << nom_frontiere << "] Boundary conditions on c and mutilde are compatible: NEUMANN on c / DIRICHLET on mutilde." << finl;
          else
            {
              Cerr << "     [" << nom_frontiere << "] Problem of compatibility between boundary conditions." << finl;
              Cerr << "         Please check your dataset and the physics behind Cahn-Hilliard equation." << finl;
              Cerr << "         Stop..." << finl;
              exit();
            }
          Cerr << "         ↪" << finl;
        }
      // 2bis) c courant imposé (pas possible sur c)
      else if (sub_type(Courant_impose,cl_c.valeur()))
        {
          Cerr << "     [" << nom_frontiere << "] Impossible to impose a current on the concentration." << finl;
          Cerr << "         Please check your dataset and the physics behind Cahn-Hilliard equation." << finl;
          Cerr << "         Stop..." << finl;
          exit();
        }
      // 3) c Dirichlet
      else if (sub_type(Scalaire_impose_paroi,cl_c.valeur()))
        {
          Cerr << "---------↧      " << finl;
          if (sub_type(Scalaire_impose_paroi,cl_mutilde.valeur()))
            {
              Cerr << "     [" << nom_frontiere << "] Boundary conditions are DIRICHLET on c / DIRICHLET on mutilde." << finl;
              Cerr << "         To keep consistance of the Cahn-Hilliard equation, BC on mutilde are computed again "
                   "depending on BC on c." << finl;
 
              Cerr << "         CANNOT CHANGE VALUES THAT ARE NOT UNIFORM. Stop..." << finl;
              exit();
 
            }
          else if (sub_type(Neumann_paroi,cl_mutilde.valeur()))
            Cerr << "     [" << nom_frontiere << "] Boundary conditions on c and mutilde are compatible: DIRICHLET on c / NEUMANN on mutilde." << finl;
          else
            {
              Cerr << "     [" << nom_frontiere << "] Problem of compatibility between boundary conditions." << finl;
              Cerr << "         Please check your dataset and the physics behind Cahn-Hilliard equation." << finl;
              Cerr << "         Stop..." << finl;
              exit();
            }
          Cerr << "         ↪" << finl;
        }
      else
        {
          Cerr << "---------↧      " << finl;
          Cerr << "     [" << nom_frontiere << "] Problem of compatibility between boundary conditions." << finl;
          Cerr << "         Please check your dataset and the physics behind Cahn-Hilliard equation." << finl;
          Cerr << "         Stop..." << finl;
          exit();
        }
    }
  return 1;
}