Fermeture_Phase_field_base.h

/****************************************************************************
* Copyright (c) 2025, 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.
*
*****************************************************************************/
 
#ifndef Fermeture_Phase_field_base_included
#define Fermeture_Phase_field_base_included
 
#include <Table.h>
 
#define PI 3.14159265
 
class Fermeture_Phase_field_base : public Objet_U
{
  Declare_base(Fermeture_Phase_field_base);
public:
  inline int get_type_systeme_naire() const { return type_systeme_naire_; }
  inline int get_type_potentiel_analytique() const { return type_potentiel_analytique_; }
  inline int get_nb_constituants() const { return nb_constituants_; }
  inline int get_type_kappa() const { return type_kappa_; }
  inline int get_type_kappa_auto_diffusion() const { return type_kappa_auto_diffusion_; }
  inline int get_kappa_ind() const { return kappa_ind_; }
 
  inline double get_mult_kappa() const { return mult_kappa_; }
  inline double get_kappa() const { return kappa_; }
  inline double get_alpha() const { return alpha_; }
  inline double get_beta() const { return beta_; }
  inline double get_molarVolume() const { return molarVolume_; }
  inline double get_prefactor() const { return prefactor_; }
  inline double get_temperature() const { return temperature_; }
 
  inline const DoubleVect& get_alphaMatrix() const { return alphaMatrix_; }
  inline const DoubleVect& get_betaMatrix() const { return betaMatrix_; }
  inline const DoubleVect& get_kappaMatrix() const { return kappaMatrix_; }
 
  double call_kappa_func_c(const double d) const { return (this->*kappa_func_c)(d); }
  DoubleTab call_kappa_func_c_naire(const DoubleTab& d) const { return (this->*kappa_func_c_naire)(d); }
 
  double (Fermeture_Phase_field_base::*kappa_func_c)(const double) const;
  DoubleTab (Fermeture_Phase_field_base::*kappa_func_c_naire)(const DoubleTab&) const;
 
  double kappa_func_c_defaut(const double) const;
  double kappa_func_c_general(const double) const;
  DoubleTab kappa_func_c_naire_general(const DoubleTab&) const;
  DoubleTab kappa_func_auto_diffusion(const DoubleTab&) const;
  DoubleTab kappa_func_c_naire_defaut(const DoubleTab&) const;
 
  double call_dWdc(const double d) const { return (this->*dWdc)(d); }
  DoubleTab call_dWdc_naire(const DoubleTab& d) const { return (this->*dWdc_naire)(d); }
 
  double (Fermeture_Phase_field_base::*dWdc)(const double) const;
  DoubleTab (Fermeture_Phase_field_base::*dWdc_naire)(const DoubleTab&) const;
 
  double dWdc_defaut(const double) const;
  double dWdc_general(const double) const;
  DoubleTab dWdc_naire_general(const DoubleTab&) const;
  DoubleTab dWdc_naire_defaut(const DoubleTab&) const;
  DoubleTab dWdc_naire_analytique_ternaire(const DoubleTab&) const;
  DoubleTab dWdc_naire_analytique_quaternaire(const DoubleTab&) const;
 
protected:
  int nb_constituants_ = -1, type_systeme_naire_ = -123;
  int type_kappa_ = -1, type_alpha_rotation_ = -123;
  int type_potentiel_analytique_ = -123, type_kappa_auto_diffusion_ = -1, kappa_ind_= -123;
 
  double alpha_ = -123., beta_ = -123., kappa_ = -123.;
  double alpha_ref_ = -123., angle_alphaMatrix_ = -123.;
  double diagonal_coeff_ = -123., mult_kappa_ = -123.;
  double temperature_ = -123., molarVolume_ = -123., prefactor_ = -123.;
 
  Table potentiel_chimique_expr_, kappa_forme_expr_;
 
  DoubleVect alphaMatrix_, kappaMatrix_, betaMatrix_;
  DoubleVect eq1_, eq2_, coeff_auto_diffusion_;
  DoubleVect angle_psi_, angle_phi_, x0Eq_, x1Eq_, x2Eq_, a0Eq_, a1Eq_, a2Eq_;
};
 
inline double Fermeture_Phase_field_base::dWdc_general(const double c) const
{
  return potentiel_chimique_expr_.val(c, 0);
}
 
inline DoubleTab Fermeture_Phase_field_base::dWdc_naire_general(const DoubleTab& c) const
{
  DoubleTab dWdc_tab(c);
  dWdc_tab = 0;
  for (int i = 0; i < c.dimension(0); i++)
    dWdc_tab(i, 0) = potentiel_chimique_expr_.val(c(i, 0), 0.);
  return dWdc_tab;
}
 
inline double Fermeture_Phase_field_base::kappa_func_c_general(const double c) const
{
  return std::max(mult_kappa_ * kappa_ * kappa_forme_expr_.val(c, 0), 0.);
}
 
inline DoubleTab Fermeture_Phase_field_base::kappa_func_c_naire_general(const DoubleTab& c) const
{
  DoubleTab kappaMatrix_tab(c.dimension(0), 1);
  for (int i = 0; i < c.dimension(0); i++)
    kappaMatrix_tab(i, 0) = std::max(mult_kappa_ * kappa_ * kappa_forme_expr_.val(c(i, 0), 0), 0.);
  return kappaMatrix_tab;
}
 
inline double Fermeture_Phase_field_base::dWdc_defaut(const double c) const
{
  return 4 * c * (c + 0.5) * (c - 0.5);
  //return 2*c*(1-c)*(1-2*c);
  //return 2*(c-0)*(c-0.68)*(2*c-0-0.68);
}
 
inline DoubleTab Fermeture_Phase_field_base::dWdc_naire_defaut(const DoubleTab& c) const
{
  DoubleTab dWdc_(c);
  dWdc_ = 0;
  for (int j = 0; j < c.line_size(); j++)
    {
      for (int i = 0; i < c.dimension(0); i++)
        {
          dWdc_(i, j) = 2 * (c(i, j) - eq1_(j)) * (c(i, j) - eq2_(j)) * (2 * c(i, j) - eq1_(j) - eq2_(j));
        }
    }
  return dWdc_;
}
 
inline DoubleTab Fermeture_Phase_field_base::dWdc_naire_analytique_ternaire(const DoubleTab& c) const
{
  DoubleTab dWdc_(c);
  dWdc_ = 0;
  DoubleTab P(c.dimension(0), 2);
  P = 0;
 
  DoubleTab dPdxU(P);
  DoubleTab dPdxZr(P);
 
  //here j defines oxide and metal
  for (int j = 0; j < 2; j++)
    {
      for (int i = 0; i < c.dimension(0); i++)
        {
          P(i, j) = pow(((cos(angle_psi_(j) * PI / 180.0) * (c(i, 0) - x0Eq_(j)) + sin(angle_psi_(j) * PI / 180.0) * (c(i, 1) - x1Eq_(j))) / (a0Eq_(j))), 2)
                    + pow(((-sin(angle_psi_(j) * PI / 180.0) * (c(i, 0) - x0Eq_(j)) + cos(angle_psi_(j) * PI / 180.0) * (c(i, 1) - x1Eq_(j))) / (a1Eq_(j))), 2);
        }
    }
 
  //here j defines oxide and metal
  for (int j = 0; j < 2; j++)
    {
      for (int i = 0; i < c.dimension(0); i++)
        {
          dPdxU(i, j) = 2 * (c(i, 0) - x0Eq_(j)) * ((pow(cos(angle_psi_(j) * PI / 180.0) / (a0Eq_(j)), 2) + pow(sin(angle_psi_(j) * PI / 180.0) / (a1Eq_(j)), 2)))
                        + 2 * (c(i, 1) - x1Eq_(j)) * ((cos(angle_psi_(j) * PI / 180.0) * sin(angle_psi_(j) * PI / 180.0)) / (pow(a0Eq_(j), 2)) - (sin(angle_psi_(j) * PI / 180.0) * cos(angle_psi_(j) * PI / 180.0)) / (pow(a1Eq_(j), 2)));
          dPdxZr(i, j) = 2 * (-(c(i, 0) - x0Eq_(j)) * sin(angle_psi_(j) * PI / 180.0) + (c(i, 1) - x1Eq_(j)) * cos(angle_psi_(j) * PI / 180.0)) * cos(angle_psi_(j) * PI / 180.0) / pow(a1Eq_(j), 2)
                         + 2 * ((c(i, 0) - x0Eq_(j)) * cos(angle_psi_(j) * PI / 180.0) + (c(i, 1) - x1Eq_(j)) * sin(angle_psi_(j) * PI / 180.0)) * sin(angle_psi_(j) * PI / 180.0) / pow(a0Eq_(j), 2);
 
        }
    }
  for (int i = 0; i < c.dimension(0); i++)
    {
      dWdc_(i, 0) = prefactor_ * (P(i, 0) * dPdxU(i, 1) + P(i, 1) * dPdxU(i, 0));
      dWdc_(i, 1) = prefactor_ * (P(i, 0) * dPdxZr(i, 1) + P(i, 1) * dPdxZr(i, 0));
    }
 
  return dWdc_;
}
 
inline DoubleTab Fermeture_Phase_field_base::dWdc_naire_analytique_quaternaire(const DoubleTab& c) const
{
  DoubleTab dWdc_(c);
  dWdc_ = 0;
  DoubleTab P(c.dimension(0), 2);
  P = 0;
  DoubleTab dPdx0(P);
  DoubleTab dPdx1(P);
  DoubleTab dPdx2(P);
 
  for (int j = 0; j < 2; j++)
    {
      for (int i = 0; i < c.dimension(0); i++)
        {
          P(i, j) = pow(
                      ((cos(angle_psi_(j) * PI / 180.0) * (c(i, 0) - x0Eq_(j)) + sin(angle_psi_(j) * PI / 180.0) * cos(angle_phi_(j) * PI / 180.0) * (c(i, 1) - x1Eq_(j))
                        + sin(angle_psi_(j) * PI / 180.0) * sin(angle_phi_(j) * PI / 180.0) * (c(i, 2) - x2Eq_(j))) / (a0Eq_(j))), 2)
                    + pow(
                      ((-sin(angle_psi_(j) * PI / 180.0) * (c(i, 0) - x0Eq_(j)) + cos(angle_psi_(j) * PI / 180.0) * cos(angle_phi_(j) * PI / 180.0) * (c(i, 1) - x1Eq_(j))
                        + cos(angle_psi_(j) * PI / 180.0) * sin(angle_phi_(j) * PI / 180.0) * (c(i, 2) - x2Eq_(j))) / (a1Eq_(j))), 2)
                    + pow(((-sin(angle_phi_(j) * PI / 180.0) * (c(i, 1) - x1Eq_(j)) + cos(angle_phi_(j) * PI / 180.0) * (c(i, 2) - x2Eq_(j))) / (a2Eq_(j))), 2);
        }
    }
  for (int j = 0; j < 2; j++)
    {
      for (int i = 0; i < c.dimension(0); i++)
        {
          dPdx0(i, j) = 2 * (c(i, 0) - x0Eq_(j)) * ((pow(cos(angle_psi_(j) * PI / 180.0) / (a0Eq_(j)), 2) + pow(sin(angle_psi_(j) * PI / 180.0) / (a1Eq_(j)), 2)))
                        + 2 * (c(i, 1) - x1Eq_(j))
                        * ((cos(angle_psi_(j) * PI / 180.0) * sin(angle_psi_(j) * PI / 180.0) * cos(angle_phi_(j) * PI / 180.0)) / (pow(a0Eq_(j), 2))
                           - (cos(angle_psi_(j) * PI / 180.0) * sin(angle_psi_(j) * PI / 180.0) * cos(angle_phi_(j) * PI / 180.0)) / (pow(a1Eq_(j), 2)))
                        + 2 * (c(i, 2) - x2Eq_(j))
                        * ((cos(angle_psi_(j) * PI / 180.0) * sin(angle_psi_(j) * PI / 180.0) * sin(angle_phi_(j) * PI / 180.0)) / (pow(a0Eq_(j), 2))
                           - (cos(angle_psi_(j) * PI / 180.0) * sin(angle_psi_(j) * PI / 180.0) * sin(angle_phi_(j) * PI / 180.0)) / (pow(a1Eq_(j), 2)));
          dPdx1(i, j) = 2 * (c(i, 0) - x0Eq_(j))
                        * ((sin(angle_psi_(j) * PI / 180.0) * cos(angle_phi_(j) * PI / 180.0) * cos(angle_phi_(j) * PI / 180.0)) / (pow(a0Eq_(j), 2))
                           - (sin(angle_psi_(j) * PI / 180.0) * cos(angle_phi_(j) * PI / 180.0) * cos(angle_phi_(j) * PI / 180.0)) / (pow(a1Eq_(j), 2)))
                        + 2 * (c(i, 1) - x1Eq_(j))
                        * ((pow(sin(angle_psi_(j) * PI / 180.0) * cos(angle_phi_(j) * PI / 180.0) / (a0Eq_(j)), 2) + pow(cos(angle_psi_(j) * PI / 180.0) * cos(angle_phi_(j) * PI / 180.0) / (a1Eq_(j)), 2))
                           + pow(sin(angle_phi_(j) * PI / 180.0) / (a2Eq_(j)), 2))
                        + 2 * (c(i, 2) - x2Eq_(j))
                        * ((pow(sin(angle_psi_(j) * PI / 180.0), 2) * cos(angle_phi_(j) * PI / 180.0) * sin(angle_phi_(j) * PI / 180.0) / (pow(a0Eq_(j), 2)))
                           + pow(cos(angle_psi_(j) * PI / 180.0), 2) * cos(angle_phi_(j) * PI / 180.0) * sin(angle_phi_(j) * PI / 180.0) / (pow(a1Eq_(j), 2))
                           + (sin(angle_phi_(j) * PI / 180.0) * cos(angle_phi_(j) * PI / 180.0) / (pow(a2Eq_(j), 2))));
          dPdx2(i, j) = 2 * (c(i, 0) - x0Eq_(j))
                        * ((sin(angle_psi_(j) * PI / 180.0) * sin(angle_phi_(j) * PI / 180.0) * cos(angle_psi_(j) * PI / 180.0)) / (pow(a0Eq_(j), 2))
                           - (sin(angle_psi_(j) * PI / 180.0) * sin(angle_phi_(j) * PI / 180.0) * cos(angle_psi_(j) * PI / 180.0)) / (pow(a1Eq_(j), 2)))
                        + 2 * (c(i, 1) - x1Eq_(j))
                        * ((pow(sin(angle_psi_(j) * PI / 180.0), 2) * sin(angle_phi_(j) * PI / 180.0) * cos(angle_phi_(j) * PI / 180.0) / (pow(a0Eq_(j), 2)))
                           + pow(cos(angle_psi_(j) * PI / 180.0), 2) * sin(angle_phi_(j) * PI / 180.0) * cos(angle_phi_(j) * PI / 180.0) / (pow(a1Eq_(j), 2))
                           + (-sin(angle_phi_(j) * PI / 180.0) * cos(angle_phi_(j) * PI / 180.0) / (pow(a2Eq_(j), 2))))
                        + 2 * (c(i, 2) - x2Eq_(j))
                        * (pow(sin(angle_psi_(j) * PI / 180.0) * sin(angle_phi_(j) * PI / 180.0) / (a0Eq_(j)), 2) + pow(cos(angle_psi_(j) * PI / 180.0) * sin(angle_phi_(j) * PI / 180.0) / (a1Eq_(j)), 2)
                           + pow(cos(angle_phi_(j) * PI / 180.0) / (a2Eq_(j)), 2));
        }
    }
  for (int i = 0; i < c.dimension(0); i++)
    {
      dWdc_(i, 0) = prefactor_ * (P(i, 0) * dPdx0(i, 1) + P(i, 1) * dPdx0(i, 0)); //+muUeq a rajouter
      dWdc_(i, 1) = prefactor_ * (P(i, 0) * dPdx1(i, 1) + P(i, 1) * dPdx1(i, 0)); //+muZrEq a rajouter
      dWdc_(i, 2) = prefactor_ * (P(i, 0) * dPdx2(i, 1) + P(i, 1) * dPdx2(i, 0)); //+muFeEq a rajouter
    }
 
  return dWdc_;
}
inline double Fermeture_Phase_field_base::kappa_func_c_defaut(const double c) const
{
  return std::max(mult_kappa_ * kappa_ * (c + 0.5) * (0.5 - c), 0.);
}
 
inline DoubleTab Fermeture_Phase_field_base::kappa_func_c_naire_defaut(const DoubleTab& c) const
{
  DoubleTab kappaMatrix_tab(c.dimension(0), 1);
  for (int i = 0; i < c.dimension(0); i++)
    kappaMatrix_tab(i, 0) = std::max(mult_kappa_ * kappa_ * (c(i, 0) + 0.5) * (0.5 - c(i, 0)), 0.);
  return kappaMatrix_tab;
 
}
 
inline DoubleTab Fermeture_Phase_field_base::kappa_func_auto_diffusion(const DoubleTab& c) const
{
  DoubleTab temp_c(c.dimension(0), c.line_size() + 1);
  DoubleTab temp_kappa(c.dimension(0), nb_constituants_, nb_constituants_); //define temp to avoid resize here
  DoubleTab kappaMat = temp_kappa;
  kappaMat = 0;
  DoubleTab kappa_ij(temp_kappa);
  DoubleTab somme_c(c.dimension(0));
  somme_c = 0;
  DoubleTab kappaMatrix_tab(c.dimension(0), nb_constituants_ * nb_constituants_);
 
  // creation d'un tableau temp_c qui contiendra toutes les concentrations des n composants non pas les n-1 seulements (n-1 eqs).
  //La derniere colonne correspondra a 1-somme sur (n-1) des concentrations
 
  for (int k = 0; k < c.line_size(); k++)
    for (int l = 0; l < c.dimension(0); l++)
      {
        temp_c(l, k) = c(l, k);
        somme_c(l) += c(l, k);
        temp_c(l, c.line_size()) = 1.0 - somme_c(l);
      }
 
  // calcul de la mobilite kappa dans un doubleTab (n1,n2,n3)
  for (int k = 0; k < nb_constituants_ + 1; k++)
    {
      for (int j = 0; j < nb_constituants_; j++)
        for (int i = 0; i < nb_constituants_; i++)
          for (int l = 0; l < c.dimension(0); l++)
            {
              //            kappa_ij(l,i,j) = (j==k) ? 1-temp_c(l,j) : -temp_c(l,j);
              //            kappa_ij(l,i,j) *= (i==k) ? 1-temp_c(l,i) : -temp_c(l,i);
              //            kappa_ij(l,i,j) *= temp_c(l,k)*coeff_auto_diff(k);
              //
              if (j == k)
                kappa_ij(l, i, j) = 1 - temp_c(l, j);
              else
                kappa_ij(l, i, j) = -temp_c(l, j);
              if (i == k)
                kappa_ij(l, i, j) *= 1 - temp_c(l, i);
              else
                kappa_ij(l, i, j) *= -temp_c(l, i);
              kappa_ij(l, i, j) *= temp_c(l, k) * coeff_auto_diffusion_(k);
            }
      kappaMat += kappa_ij;
    }
 
  //retour vers un DoubleTab (n1,n2)
  for (int j = 0; j < nb_constituants_; j++)
    for (int i = 0; i < nb_constituants_; i++)
      for (int l = 0; l < c.dimension(0); l++)
        kappaMatrix_tab(l, i + nb_constituants_ * j) = kappaMat(l, i, j);
 
  kappaMatrix_tab *= molarVolume_ / prefactor_ * 8.314 * temperature_;
  return kappaMatrix_tab;
 
}
 
#endif /* Fermeture_Phase_field_base_included */