Domaine_ALE.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 <Scatter.h>
#include <MD_Vector_tools.h>
#include <MD_Vector_std.h>
#include <Debog.h>
#include <Domaine_ALE.h>
#include <Probleme_base.h>
#include <Schema_Temps_base.h>
#include <Domaine_VEF.h>
#include <Domaine_EF.h>
#include <Domaine_Poly_base.h>
#include <Motcle.h>
#include <EFichier.h>
#include <Champ_front_ALE.h>
#include <Ch_front_input_ALE.h>
#include <Champ_front_ALE_Beam.h>
#include <Equation_base.h>
#include <Noms.h>
#include <Navier_Stokes_std.h>
#include <Operateur_Diff.h>
#include <Operateur_Grad.h>
#include <communications.h>
#include <Faces.h>
#include <CL_Types_include.h>
#include <Entree_fluide_vitesse_imposee_ALE.h>
#include <Interprete_bloc.h>
#include <NettoieNoeuds.h>
#include <TRUST_2_MED.h>
#include <Matrix_tools.h>
#include <Array_tools.h>
#include <Champ_front_ALE_lag.h>
#include <EcritureLectureSpecial.h>
#include <Avanc.h>
#include <TRUST_2_PDI.h>
 
#include <medcoupling++.h>
#ifdef MEDCOUPLING_
#include <MEDCouplingMemArray.hxx>
using namespace MEDCoupling;
#endif
 
 
 
 
// (ALE convection grammar is declared as `convection_ale` near
// Op_Conv_ALE_VEF::Implemente_instanciable — the earlier
// `bloc_convection_ALE convection_deriv ALE NO_BRACE` + recursive
// `operateur convection_deriv` shape modelled `convection { ALE
// muscl }` but tripped on `convection { ALE muscl vanleer 2 }`
// because the inner muscl had no attrs to absorb the limiter and
// order. Replaced by the flat positional shape used by
// convection_generic.)
 
 
// Large penalization coefficient to enforce Dirichlet-like constraints in the Laplacian system
namespace { static constexpr double ALE_PENALIZATION = 1.e14; }
 
Implemente_instanciable_sans_constructeur(Domaine_ALE, "Domaine_ALE", Domaine);
// XD domaine_ale domaine domaine_ale INHERITS_BRACE Domain with nodes at the interior of the domain which are
// displaced in an arbitrarily prescribed way thanks to ALE (Arbitrary Lagrangian-Eulerian) description.
// NL2 Keyword to specify that the domain is mobile following the displacement of some of its boundaries.
 
 
Domaine_ALE::Domaine_ALE()
{
  clear();
}
 
void Domaine_ALE::clear()
{
  Domaine::clear();
  deformable_ = true;
 
  dt_ = 0.;
  ALE_mesh_velocity.reset();
  vf.reset();
  som_faces_bords.reset();
  solv = SolveurSys();
  les_champs_front = TRUST_Vector<OWN_PTR(Champ_front_base)>();
  nb_bords_ALE = 0;
  les_bords_ALE.vide();
  resumption = 0;
  eq = OBS_PTR(Equation_base)();
  les_elems_extrait_surf_reference_.reset();
  extrait_surf_dom_deformable_ = false;
  str_mesh_model = Structural_dynamic_mesh_model();
  meshMotionModel_ = 0;
  UpdateTheGrid = true;
  Coupling_ICoCo_method = -1;
  volumes_old_.reset();
  volumes_entrelaces_old_.reset();
  coord_old_.reset();
 
  beam_coupling_.clear();
  projection_manager_.clear();
}
 
Sortie& Domaine_ALE::printOn(Sortie& os) const
{
  Domaine::printOn(os);
  return os;
}
 
Entree& Domaine_ALE::readOn(Entree& is)
{
  Domaine::readOn(is);
  return is;
}
 
 
void Domaine_ALE::mettre_a_jour(double temps, Domaine_dis_base& le_domaine_dis,
                                Probleme_base& pb)
{
  invalide_octree();
 
  // Initialization of the ICoCo coupling method on first call
  if (Coupling_ICoCo_method == -1)
    {
      int res = getCouplingMethod();
      Cout << "Initialization of Coupling_ICoCo_Method " << res << finl;
    }
 
  //--------------------------------------------------------------------------
  // ICoCo implicit sub-iteration handling: save/reset state at iteration boundaries
  //--------------------------------------------------------------------------
  int iCoCoImplicitIteration = -1; // -1 = no ICoCo implicit coupling
  bool resetCoordinates = false;
 
  if ((eq.valeur()).probleme().checkOutputIntEntry("iCoCoImplicitIteration"))
    {
      iCoCoImplicitIteration = (eq.valeur()).probleme().getOutputIntValue("iCoCoImplicitIteration");
      Cerr << "Domaine_ALE::mettre_a_jour, iCoCoImplicitIteration: "
           << iCoCoImplicitIteration << finl;
      if (meshMotionModel_ == 1)
        str_mesh_model.iCoCoImplicitIteration = iCoCoImplicitIteration;
    }
 
  if (meshMotionModel_ == 1)
    {
      str_mesh_model.acceleratedSolution = 0; // Reset to default each step
      if ((eq.valeur()).probleme().checkOutputIntEntry("acceleratedGridMode"))
        str_mesh_model.acceleratedSolution =
          (eq.valeur()).probleme().getOutputIntValue("acceleratedGridMode");
    }
 
  if (iCoCoImplicitIteration == 0)
    {
      saveSommetsCoordinates(); // Save coordinates at first ICoCo sub-iteration
      if (meshMotionModel_ == 1) str_mesh_model.saveBeginning_of_stepVariables();
    }
  else if (iCoCoImplicitIteration > 0)
    {
      resetSommetsCoordinates(); // Reset to beginning-of-step for subsequent iterations
      resetCoordinates = true;
      if (meshMotionModel_ == 1)
        {
          if (str_mesh_model.acceleratedSolution == 0)
            str_mesh_model.resetVariablestoBeginning_of_step();
          else
            str_mesh_model.gridTime = str_mesh_model.gridTime_n; // Accelerated mode: reset time only
        }
    }
 
  compute_vertex_velocities(temps, pb);
 
  if (!(mesh_update_required_ || resetCoordinates)) return;
 
  Cout << "Updating the mesh!" << finl;
  const double dt = pb.schema_temps().pas_de_temps();
  update_vertex_coordinates(dt);
 
  // Recompute all geometric metrics after mesh motion
  updateMetrics(le_domaine_dis, pb);
 
  // Recompute face velocities from updated vertex velocities
  Domaine_VF& le_dom_VF = ref_cast(Domaine_VF, le_domaine_dis);
  calculer_vitesse_faces(ALE_mesh_velocity, le_dom_VF.nb_faces(),
                         le_dom_VF.nb_som_face(), le_dom_VF.face_sommets());
 
  update_post_fields(temps, dt);
 
  // Enforce exact periodicity at periodic face pairs.
  // Calling verifie_valeurs_cl() here (before the NS operators
  // run) keeps the periodic constraint tight at every time step.
  if (eq)
    {
      const Domaine_Cl_dis_base& zcl = eq.valeur().domaine_Cl_dis();
      bool has_periodic = false;
      for (int cl = 0; cl < zcl.nb_cond_lim() && !has_periodic; cl++)
        if (sub_type(Periodique, zcl.les_conditions_limites(cl).valeur()))
          has_periodic = true;
      if (has_periodic)
        eq.valeur().inconnue().verifie_valeurs_cl();
    }
}
 
 
void Domaine_ALE::updateMetrics(Domaine_dis_base& le_domaine_dis, Probleme_base& pb)
{
  if (!eq) eq = pb.equation(0);
  metrics_updater_.update(le_domaine_dis, pb, *this);
}
 
 
void Domaine_ALE::update_after_post(double temps)
{
  // Compute and write modal fluid forces for all projection boundaries
  projection_manager_.update(temps, nb_bords_ALE, les_bords_ALE, nom_du_cas(), getEquation());
}
 
 
 
void Domaine_ALE::update_ALE_projection(double temps, Nom& name_ALE_boundary_projection,
                                        Champ_front_ALE_projection& field_ALE_projection,
                                        int nb_mode)
{
  projection_manager_.update(temps, name_ALE_boundary_projection, field_ALE_projection,
                             nb_mode, nb_bords_ALE, les_bords_ALE,
                             nom_du_cas(), getEquation());
}
 
void Domaine_ALE::update_ALE_projection(const double temps)
{
  projection_manager_.update(temps, nb_bords_ALE, les_bords_ALE, nom_du_cas(), getEquation());
}
 
 
 
void Domaine_ALE::update_post_fields(double temps, double dt)
{
  const DoubleTab& ALEMeshVelocity = vitesse_faces();
 
  if (ALEMeshVelocity_)
    {
      ALEMeshVelocity_->valeurs() = ALEMeshVelocity;
      ALEMeshVelocity_->mettre_a_jour(temps);
    }
 
  if (ALEMeshTotalDisplacement_)
    {
      // Integrate mesh velocity: disp_{n+1} = disp_n + dt * v_mesh
      const int D = dimension;
      const int N = ALEMeshVelocity.size() / D;
      DoubleTab& disp = ALEMeshTotalDisplacement_->valeurs();
      for (int i = 0; i < N; i++)
        for (int d = 0; d < D; d++)
          disp(i, d) += ALEMeshVelocity(i, d) * dt;
      disp.echange_espace_virtuel();
      ALEMeshTotalDisplacement_->mettre_a_jour(temps);
    }
 
  if (meshMotionModel_ == 1)
    {
      if (ALEMeshStructuralPressure_)
        ALEMeshStructuralPressure_->valeurs() = getMeshPbPressure();
      if (ALEMeshStructuralVonMises_)
        ALEMeshStructuralVonMises_->valeurs() = getMeshPbVonMises();
      if (ALEMeshStructuralForces_)
        ALEMeshStructuralForces_->valeurs() = getMeshPbForceFace();
    }
}
 
 
void Domaine_ALE::validate_ale_setup(const Probleme_base& pb) const
{
  // Ensure Neumann boundaries are not also declared as moving boundaries
  const Noms& neumann = projection_manager_.neumann_boundary_names();
  for (int i = 0; i < neumann.size(); i++)
    for (int j = 0; j < nb_bords_ALE; j++)
      if (les_bords_ALE(j).le_nom() == neumann[i])
        {
          Cerr << "In 'ALE_Neumann_BC_for_grid_problem', boundary " << neumann[i]
               << " is already declared as a moving boundary in 'Imposer_vit_bords_ALE'." << finl;
          Process::exit();
        }
 
  // Check discretization compatibility
  if (!(sub_type(Domaine_VEF, pb.domaine_dis()) ||
        pb.discretisation().is_PolyMAC_MPFA() ||
        pb.discretisation().is_ef()))
    {
      Cerr << "Discretization not recognized by ALE!" << finl;
      Cerr << "Use VEFPreP1B, PolyMAC_P0 or EF and restart." << finl;
      Process::exit();
    }
 
  // Check that moving boundaries use the correct boundary condition type
  const Domaine_Cl_dis_base& dom_cl = pb.equation(0).domaine_Cl_dis();
  for (int j = 0; j < nb_bords_ALE; j++)
    {
      const Nom& le_nom = les_bords_ALE(j).le_nom();
      const int rang    = rang_frontiere(le_nom);
      if (rang >= faces_bord().size()) continue; // internal boundary, skip
      const Cond_lim& la_cl = dom_cl.les_conditions_limites(rang);
      if (!sub_type(Entree_fluide_vitesse_imposee_ALE, la_cl.valeur()))
        Cerr << "Mobile ALE boundary: replace " << la_cl->que_suis_je()
             << " on " << le_nom << " with Frontiere_ouverte_vitesse_imposee_ALE" << finl;
    }
 
  if (les_champs_front.size() == 0)
    Cerr << "Warning: ALE simulation without any moving boundary!" << finl;
}
 
void Domaine_ALE::initialiser(double temps, Domaine_dis_base& le_domaine_dis,
                              Probleme_base& pb)
{
  const bool initialized = (vf.size_array() != 0);
  if (initialized) return;
 
  Cerr << "Domaine_ALE::initialize" << finl;
  validate_ale_setup(pb);
  invalide_octree();
  eq = pb.equation(0);
 
  Domaine_VF& le_dom_VF = ref_cast(Domaine_VF, le_domaine_dis);
  const int nb_faces     = le_dom_VF.nb_faces();
  const int nb_faces_tot = le_dom_VF.nb_faces_tot();
  const int nb_som_face  = le_dom_VF.nb_som_face();
  IntTab& face_sommets   = le_dom_VF.face_sommets();
 
  // Initialize vertex and face velocity arrays with proper parallel distribution
  const MD_Vector& mds = md_vector_sommets();
  MD_Vector_tools::creer_tableau_distribue(mds, laplacian_rhs_);
  ALE_mesh_velocity.resize(0, dimension);
  MD_Vector_tools::creer_tableau_distribue(mds, ALE_mesh_velocity);
  vf.resize(nb_faces, dimension);
  MD_Vector_tools::creer_tableau_distribue(le_dom_VF.md_vector_faces(), vf);
 
  // Snapshot initial metrics for the first time step's volume corrections
  coord_old_              = sommets_;
  volumes_old_            = le_dom_VF.volumes();
  volumes_entrelaces_old_ = le_dom_VF.volumes_entrelaces();
  if (pb.discretisation().is_ef())
    volumes_sommets_thilde_old_ = ref_cast(Domaine_EF, le_dom_VF).volumes_sommets_thilde();
 
  if (meshMotionModel_ == 1)
    {
      // Initialize the fictitious structural problem for mesh motion
      str_mesh_model.initMfrontBehaviour();
      str_mesh_model.initDynamicMeshProblem(
        temps, nb_som(), nb_elem(), nb_faces,
        mds, md_vector_elements(), le_dom_VF.md_vector_faces());
    }
 
  if (lagrangian_velocity_imposed_)
    {
#ifdef MEDCOUPLING_
      const Probleme_base& pb_v = ref_cast(Probleme_base, interprete().objet(external_pb_name_));
      const Champ_base& champ   = pb_v.get_champ(external_velocity_field_name_);
      build_lagrangian_vertex_mapping(pb_v, champ);
#else
      Cerr << "Domaine_ALE::initialiser : MEDCoupling required for imposed Lagrangian velocity." << finl;
      Process::exit();
#endif
    }
 
  // Compute initial mesh velocities (using resumption time if restarting)
  if (resumption)
    compute_vertex_velocities(pb.schema_temps().temps_courant(), pb);
  else
    compute_vertex_velocities(temps, pb);
 
  calculer_vitesse_faces(ALE_mesh_velocity, nb_faces_tot, nb_som_face, face_sommets);
 
  // Initialize Ch_front_input_ALE fields if present
  for (int i = 0; i < les_champs_front.size(); i++)
    {
      Champ_front_base& field = les_champs_front[i].valeur();
      if (sub_type(Ch_front_input_ALE, field))
        {
          const Champ_Inc_base& vit = ref_cast(Champ_Inc_base, pb.get_champ("VITESSE"));
          ref_cast(Ch_front_input_ALE, field).initialiser(temps, vit);
        }
    }
}
 
DoubleTab& Domaine_ALE::calculer_vitesse_faces(DoubleTab& vit_maillage, int nb_faces,
                                               int nb_som_face, IntTab& face_sommets)
{
  vf = 0.;
  // If the mesh velocity is essentially zero, face velocity stays zero (default).
  if (vit_maillage.mp_max_abs_vect() >= 1.e-12)
    for (int j = 0; j < dimension; j++)
      for (int i = 0; i < nb_faces; i++)
        {
          for (int s = 0; s < nb_som_face; s++)
            vf(i, j) += vit_maillage(face_sommets(i, s), j);
          vf(i, j) /= nb_som_face;
        }
 
  vf.echange_espace_virtuel();
  Debog::verifier("Domaine_ALE::calculer_vitesse_faces", vf);
  return vf;
}
 
void Domaine_ALE::set_dt(double& dt) { dt_ = dt; }
 
 
void Domaine_ALE::validateTimeStep()
{
  Domaine::validateTimeStep();
  if (eq)
    {
      // Snapshot metrics so that volume corrections use V_n while the solver uses V_{n+1}
      Domaine_VF& le_dom_VF = ref_cast(Domaine_VF, eq->domaine_dis());
      coord_old_ = sommets_;
      volumes_old_ = le_dom_VF.volumes();
      if (!eq->probleme().discretisation().is_ef())
        volumes_entrelaces_old_ = le_dom_VF.volumes_entrelaces();
      else
        volumes_sommets_thilde_old_ = ref_cast(Domaine_EF, le_dom_VF).volumes_sommets_thilde();
    }
  resumption = 0;
}
 
 
void Domaine_ALE::ajouter_correctif_volumique(const DoubleTab& v, const DoubleTab& derivee,
                                              double dt, DoubleTab& deriveeALE) const
{
  // ALE volume correction: (V_n / V_{n+1}) * (U_n / dt) + derivee
  const DoubleVect& volumes_entrelaces =
    ref_cast(Domaine_VF, eq->domaine_dis()).volumes_entrelaces();
  const int nf_tot = v.dimension_tot(0);
  for (int f = 0; f < nf_tot; f++)
    for (int d = 0; d < dimension; d++)
      deriveeALE(f, d) = v(f, d) * volumes_entrelaces_old_(f) / volumes_entrelaces(f) / dt
                         + derivee(f, d);
 
  // Enforce periodicity: the V_old/V_new ratio can differ slightly between periodic
  // partner faces when the ALE mesh deforms non-uniformly, breaking the symmetry that
  // Op_Div_VEFP1B_Elem::ajouter requires. Average each periodic pair so that both
  // partners carry the same value, consistent with how verifie_valeurs_cl treats velocity.
  if (eq)
    {
      const Domaine_Cl_dis_base& zcl = eq.valeur().domaine_Cl_dis();
      for (int i = 0; i < zcl.nb_cond_lim(); i++)
        {
          const Cond_lim_base& la_cl = zcl.les_conditions_limites(i).valeur();
          if (sub_type(Periodique, la_cl))
            {
              const Periodique& la_cl_perio = ref_cast(Periodique, la_cl);
              const Front_VF& le_bord = ref_cast(Front_VF, la_cl.frontiere_dis());
              int ndeb = le_bord.num_premiere_face();
              int nfin = ndeb + le_bord.nb_faces();
              for (int f = ndeb; f < nfin; f++)
                {
                  int fass = la_cl_perio.face_associee(f - ndeb) + ndeb;
                  for (int d = 0; d < dimension; d++)
                    {
                      double moy = 0.5 * (deriveeALE(f, d) + deriveeALE(fass, d));
                      deriveeALE(f, d) = moy;
                      deriveeALE(fass, d) = moy;
                    }
                }
            }
        }
    }
}
 
void Domaine_ALE::apply_old_to_new_volume_scaling(DoubleTab& inco,
                                                  const Domaine_dis_base& ddb) const
{
  const Domaine_VF& dvf = ref_cast(Domaine_VF, ddb);
  const int I = inco.dimension_tot(0);
  const int N = inco.dimension(1);
 
  if (I == nb_elem_tot())
    {
      // Cell-centered: scale by V_n / V_{n+1}
      const DoubleVect& volumes = dvf.volumes();
      for (int e = 0; e < I; e++)
        for (int n = 0; n < N; n++)
          inco(e, n) *= volumes_old_(e) / volumes(e);
    }
  else if (I == dvf.nb_faces_tot())
    {
      const DoubleVect& ve = dvf.volumes_entrelaces();
      for (int f = 0; f < I; f++)
        for (int n = 0; n < N; n++)
          inco(f, n) *= volumes_entrelaces_old_(f) / ve(f);
    }
  else if (I == dvf.nb_som_tot())
    {
      const DoubleVect& vs = ref_cast(Domaine_EF, dvf).volumes_sommets_thilde();
      for (int s = 0; s < I; s++)
        for (int n = 0; n < N; n++)
          inco(s, n) *= volumes_sommets_thilde_old_(s) / vs(s);
    }
  else
    {
      Cerr << "Domaine_ALE::apply_old_to_new_volume_scaling: incompatible dimension "
           << I << " (expected " << nb_elem_tot() << ", " << dvf.nb_faces_tot() << ")." << finl;
      Process::exit();
    }
}
 
//////////////////////////////////////////////////////////////////////////////
// Structural dynamics field accessors (forwarded to str_mesh_model)
//////////////////////////////////////////////////////////////////////////////
 
void Domaine_ALE::resumptionStructuralDynamicsMesh(
  double tinit, DoubleTab& u_n, DoubleTab& v_n, DoubleTab& a_n,
  DoubleTab& x_n, DoubleTab& B0_n, DoubleTab& Ft_n, DoubleTab& Stress_n)
{
  str_mesh_model.resumptionMesh(tinit, u_n, v_n, a_n, x_n, B0_n, Ft_n, Stress_n);
}
 
const DoubleVect& Domaine_ALE::getMeshPbPressure()              const { return str_mesh_model.getMeshPbPressure(); }
const DoubleVect& Domaine_ALE::getMeshPbVonMises()              const { return str_mesh_model.getMeshPbVonMises(); }
const DoubleTab&  Domaine_ALE::getMeshPbForceFace()             const { return str_mesh_model.getMeshPbForceFace(); }
const DoubleTab&  Domaine_ALE::getMeshDisplacement()            const { return str_mesh_model.getMeshDisplacement(); }
const DoubleTab&  Domaine_ALE::getMeshVelocity()                const { return str_mesh_model.getMeshVelocity(); }
const DoubleTab&  Domaine_ALE::getMeshAcceleration()            const { return str_mesh_model.getMeshAcceleration(); }
const DoubleTab&  Domaine_ALE::getMeshPosition()                const { return str_mesh_model.getMeshPosition(); }
const DoubleTab&  Domaine_ALE::getMeshReferenceConfiguration()  const { return str_mesh_model.getMeshReferenceConfiguration(); }
const DoubleTab&  Domaine_ALE::getMeshTransformationGradient()  const { return str_mesh_model.getMeshTransformationGradient(); }
const DoubleTab&  Domaine_ALE::getMeshStress()                  const { return str_mesh_model.getMeshStress(); }
const int&        Domaine_ALE::getMeshReferenceConfigurationNbComp() const { return str_mesh_model.getMeshReferenceConfigurationNbComp(); }
const int&        Domaine_ALE::getMeshTransformationGradientNbComp() const { return str_mesh_model.getMeshTransformationGradientNbComp(); }
const int&        Domaine_ALE::getMeshStressNbComp()            const { return str_mesh_model.getMeshStressNbComp(); }
 
//////////////////////////////////////////////////////////////////////////////
// Beam accessors  (forwarded to beam_coupling_)
//////////////////////////////////////////////////////////////////////////////
 
const Nom&       Domaine_ALE::getBeamName       (const int& i)              const { return beam_coupling_.getName(i); }
const int&       Domaine_ALE::getBeamNbModes    (const int& i)              const { return beam_coupling_.getNbModes(i); }
const int&       Domaine_ALE::getBeamNbPlanes   (const int& i)              const { return beam_coupling_.getNbPlanes(i); }
const int&       Domaine_ALE::getBeamNbBeam     ()                          const { return beam_coupling_.getNbBeam(); }
const int&       Domaine_ALE::getBeamLongitudinalAxis(const int& i)         const { return beam_coupling_.getLongitudinalAxis(i); }
const int&       Domaine_ALE::getBeamBendingDirection(const int& i, const int& idx) const { return beam_coupling_.getBendingDirection(i, idx); }
const DoubleTab& Domaine_ALE::getBeamDisplacement(const int& i, const int& j) const { return beam_coupling_.getDisplacement(i, j); }
const DoubleTab& Domaine_ALE::getBeamRotation    (const int& i, const int& j) const { return beam_coupling_.getRotation(i, j); }
 
DoubleTab& Domaine_ALE::getBeamVelocity(const int& i, const double& tps, const double& dt)
{
  // eq may be null during BC initialization at t=0 (before Domaine_ALE::initialiser runs).
  // getVelocity only dereferences eqn when dt != 0, so passing nullptr is safe here.
  Equation_base* eqn = eq ? &eq.valeur() : nullptr;
  return beam_coupling_.getVelocity(i, tps, dt, eqn, is_resumption(), nb_bords_ALE, les_bords_ALE);
}
 
void Domaine_ALE::computeFluidForceOnBeam(const int& i)
{
  beam_coupling_.computeFluidForce(i, getEquation(), is_resumption(), nb_bords_ALE, les_bords_ALE);
}
 
double Domaine_ALE::computeDtBeam(Domaine_dis_base& le_domaine_dis, const int& i)
{
  return beam_coupling_.computeDt(le_domaine_dis, i);
}
 
DoubleTab Domaine_ALE::interpolationOnThe3DSurface(
  const int& i, const double& x, const double& y, const double& z,
  const DoubleTab& u, const DoubleTab& R) const
{
  return beam_coupling_.interpolationOnThe3DSurface(i, x, y, z, u, R);
}
 
double Domaine_ALE::interpolationPhiOnThe3DSurface(
  const int& i, const double& x, const double& y, const double& z,
  const int& comp, const DoubleTab& u) const
{
  return beam_coupling_.interpolationPhiOnThe3DSurface(i, x, y, z, comp, u);
}
 
Equation_base& Domaine_ALE::getEquation()
{
  return eq;
}
 
const Equation_base& Domaine_ALE::getEquation() const
{
  return eq;
}
 
//////////////////////////////////////////////////////////////////////////////
// Checkpoint delegates
//////////////////////////////////////////////////////////////////////////////
 
std::vector<YAML_data> Domaine_ALE::data_a_sauvegarder(const Probleme_base& pb) const
{
  return checkpoint_manager_.data_a_sauvegarder(pb, *this);
}
 
int Domaine_ALE::save_additional_state(Sortie& os, const Probleme_base& pb) const
{
  if (!eq) eq = pb.equation(0);
  return checkpoint_manager_.save(os, pb, *this);
}
 
int Domaine_ALE::restore_additional_state(Entree& is, Probleme_base& pb)
{
  resumption = 1;
  if (!eq) eq = pb.equation(0);
  return checkpoint_manager_.restore(is, pb, *this);
}
 
Nom Domaine_ALE::generate_field_tag(const Champ_Inc_base& champ,
                                    const Probleme_base& pb) const
{
  return checkpoint_manager_.generate_field_tag(champ, pb);
}
 
YAML_data Domaine_ALE::make_yaml(const std::string& suffix, int nb_dim,
                                 const Probleme_base& pb) const
{
  return checkpoint_manager_.make_yaml(suffix, nb_dim, pb);
}
 
void Domaine_ALE::create_field(OWN_PTR(Champ_Inc_base)& champ, const std::string& name,
                               int nbComp, const Motcle& directive,
                               const Probleme_base& pb) const
{
 
  checkpoint_manager_.create_field(champ, name, nbComp, directive, pb,
                                   const_cast<Domaine_ALE&>(*this));
}
 
//////////////////////////////////////////////////////////////////////////////
// Resumption Coordinates
//////////////////////////////////////////////////////////////////////////////
 
 
void Domaine_ALE::resumptionCoords(DoubleTab& ValueOf_MeshCoords)
{
  coord_old_ = sommets_;
  const int nbValues = ValueOf_MeshCoords.size() / dimension;
  for (int i = 0; i < nbValues; i++)
    for (int k = 0; k < dimension; k++)
      coord(i, k) = ValueOf_MeshCoords(i, k);
}
 
 
void Domaine_ALE::reading_vit_bords_ALE(Entree& is)
{
  Motcle accolade_ouverte("{"), accolade_fermee("}"), motlu;
  Nom nomlu;
 
  is >> motlu;
  if (motlu != accolade_ouverte)
    {
      Cerr << "Error reading moving boundary velocities: expected "
           << accolade_ouverte << ", got " << motlu << finl;
      Process::exit();
    }
 
  is >> nb_bords_ALE;
  Cerr << "Number of ALE boundaries : " << nb_bords_ALE << finl;
  les_champs_front.dimensionner(nb_bords_ALE);
 
  int compteur = 0;
  while (1)
    {
      is >> nomlu;
      motlu = nomlu;
      if (motlu == accolade_fermee) break;
      Cerr << "Reading velocity imposed on " << nomlu << "..." << finl;
      les_bords_ALE.add(faces_bord()(rang_frontiere(nomlu)));
      is >> les_champs_front[compteur++];
    }
}
 
void Domaine_ALE::reading_solver_moving_mesh_ALE(Entree& is)
{
  Motcle accolade_ouverte("{"), accolade_fermee("}"), motlu;
  Nom nomlu;
 
  is >> motlu;
  if (motlu != accolade_ouverte)
    {
      Cerr << "Error reading ALE mesh solver: expected "
           << accolade_ouverte << ", got " << motlu << finl;
      Process::exit();
    }
 
  is >> solv;
  solv.nommer("ALE_solver");
 
  while (1)
    {
      is >> nomlu;
      motlu = nomlu;
      if (motlu == accolade_fermee) break;
    }
}
 
void Domaine_ALE::reading_beam_model(Entree& is)
{
  beam_coupling_.read(is, nb_bords_ALE, les_bords_ALE);
}
 
void Domaine_ALE::reading_projection_ALE_boundary(Entree& is)
{
  projection_manager_.read_projection_boundary(is);
}
 
void Domaine_ALE::reading_ALE_Neumann_BC_for_grid_problem(Entree& is)
{
  projection_manager_.read_neumann_bc(is);
}
 
void Domaine_ALE::reading_structural_dynamic_mesh_model(Entree& is)
{
  meshMotionModel_ = 1;
 
  Motcle accolade_ouverte("{"), accolade_fermee("}"), motlu;
  Nom nomlu;
  double var_double;
  std::vector<double> val_prop;
 
  is >> motlu;
  if (motlu != accolade_ouverte)
    {
      Cerr << "Error reading structural dynamic mesh model: expected "
           << accolade_ouverte << ", got " << motlu << finl;
      Process::exit();
    }
 
  while (1)
    {
      is >> nomlu;
      motlu = nomlu;
 
      if (motlu == "Mfront_library")
        { is >> nomlu; str_mesh_model.setMfrontLibraryPath(nomlu.getString()); continue; }
      if (motlu == "Mfront_model_name")
        { is >> nomlu; str_mesh_model.setMfrontModelName(nomlu.getString()); continue; }
      if (motlu == "Mfront_hypothesis")
        { is >> nomlu; str_mesh_model.setMfrontHypothesis(nomlu.getString()); continue; }
      if (motlu == "Mfront_material_property")
        {
          is >> motlu;
          if (motlu != accolade_ouverte)
            { Cerr << "Error: expected { in Mfront_material_property" << finl; Process::exit(); }
          is >> nomlu;
          std::string nom_prop(nomlu);
          val_prop.clear();
          is >> var_double;
          while (1)
            {
              val_prop.push_back(var_double);
              is >> motlu;
              if (motlu == accolade_fermee) break;
              var_double = std::stod(motlu.getString());
            }
          str_mesh_model.addMaterialProperty(nom_prop, val_prop);
          continue;
        }
      if (motlu == "Density")
        { is >> var_double; str_mesh_model.setDensity(var_double); continue; }
      if (motlu == "Inertial_Damping")
        { is >> var_double; str_mesh_model.setInertialDamping(var_double); continue; }
      if (motlu == "Time_Step_Safety_Coefficient")
        { is >> var_double; str_mesh_model.setDtSafetyCoefficient(var_double); continue; }
      if (motlu == "Grid_Dt_Min")
        { is >> var_double; str_mesh_model.setGridDtMin(var_double); continue; }
      if (motlu == "Configuration_Reset_Dt")
        { is >> var_double; str_mesh_model.setConfigurationResetDt(var_double); continue; }
      if (motlu == "Max_Added_Mass_Ratio")
        { is >> var_double; str_mesh_model.setMaxAddedMassRatio(var_double); continue; }
      if (motlu == accolade_fermee) break;
    }
 
  if (str_mesh_model.getDensity() == 0.)
    { Cerr << "Error: density not provided in structural dynamic mesh model." << finl; Process::exit(); }
}
 
//////////////////////////////////////////////////////////////////////////////
// ICoCo coupling
//////////////////////////////////////////////////////////////////////////////
 
int Domaine_ALE::getCouplingMethod()
{
  if (Coupling_ICoCo_method == -1) // First call: initialize from ICoCo or use default
    {
      Probleme_base& pb_base = (eq.valeur()).probleme();
      if (pb_base.checkOutputIntEntry("CouplingMethod"))
        {
          Coupling_ICoCo_method = pb_base.getOutputIntValue("CouplingMethod");
          Cerr << "Coupling Method set to " << Coupling_ICoCo_method << " from ICoCo" << finl;
        }
      else
        {
          Coupling_ICoCo_method = 0;
          Cerr << "Coupling Method set to " << Coupling_ICoCo_method << " (default)" << finl;
        }
    }
  return Coupling_ICoCo_method;
}
 
bool Domaine_ALE::getCouplingMethod() const
{
  return (Coupling_ICoCo_method == 1);
}
 
bool Domaine_ALE::getUpdateTheGrid()
{
  Probleme_base& pb_base = (eq.valeur()).probleme();
  if (pb_base.checkOutputIntEntry("UpdateTheGrid"))
    if (pb_base.getOutputIntValue("UpdateTheGrid") == 1)
      {
        Cerr << "Update of the fluid grid forced from ICoCo" << finl;
        UpdateTheGrid = true;
      }
  return UpdateTheGrid;
}
 
void Domaine_ALE::setUpdateTheGrid(bool val) { UpdateTheGrid = val; }
 
//////////////////////////////////////////////////////////////////////////////
// Lagrangian external velocity coupling
//////////////////////////////////////////////////////////////////////////////
 
void Domaine_ALE::associer_vitesse_lagrangienne(const Nom& nom_pb, const Nom& nom_champ)
{
  external_pb_name_            = nom_pb;
  external_velocity_field_name_ = nom_champ;
  lagrangian_velocity_imposed_  = true;
}
 
void Domaine_ALE::check_external_vertex_velocity_compatibility(const DoubleTab& valeurs) const
{
  if (valeurs.line_size() != dimension)
    {
      Cerr << "Domaine_ALE: Lagrangian velocity field must have " << dimension
           << " components per node (line_size=" << valeurs.line_size() << ")." << finl;
      Process::exit();
    }
}
 
void Domaine_ALE::build_lagrangian_vertex_mapping(const Probleme_base& pb_v,
                                                  const Champ_base& champ)
{
#ifndef MEDCOUPLING_
  Cerr << "Domaine_ALE::build_lagrangian_vertex_mapping: MEDCoupling required." << finl;
  Process::exit();
#else
  const DoubleTab& ext_coords = pb_v.domaine().les_sommets();
  const DoubleTab& ale_coords = sommets_;
  const DoubleTab& valeurs    = champ.valeurs();
 
  check_external_vertex_velocity_compatibility(valeurs);
  if (valeurs.dimension(0) != ext_coords.dimension(0))
    {
      Cerr << "Domaine_ALE: external velocity field has " << valeurs.dimension(0)
           << " nodes but domain has " << ext_coords.dimension(0) << " nodes." << finl;
      Process::exit();
    }
 
  const int nb_ext = ext_coords.dimension(0);
  const int nb_ale = ale_coords.dimension(0);
  const int dim    = dimension;
  const double tol = Objet_U::precision_geom * 100.;
 
  external_vertex_id_for_ale_vertex_.resize(nb_ale);
  external_vertex_id_for_ale_vertex_ = -1;
 
  int unmatched = 0;
  if (nb_ext == 0)
    {
      unmatched = nb_ale;
    }
  else
    {
      // Use MEDCoupling nearest-neighbor search for efficiency
      MCAuto<DataArrayDouble> ext_pts(DataArrayDouble::New());
      MCAuto<DataArrayDouble> ale_pts(DataArrayDouble::New());
      ext_pts->alloc(nb_ext, dim);
      ale_pts->alloc(nb_ale, dim);
 
      double* ext_ptr = ext_pts->getPointer();
      for (int i = 0; i < nb_ext; i++)
        for (int d = 0; d < dim; d++)
          ext_ptr[i * dim + d] = ext_coords(i, d);
 
      double* ale_ptr = ale_pts->getPointer();
      for (int i = 0; i < nb_ale; i++)
        for (int d = 0; d < dim; d++)
          ale_ptr[i * dim + d] = ale_coords(i, d);
 
      MCAuto<DataArrayIdType> closest(ext_pts->findClosestTupleId(ale_pts));
      const double tol2 = tol * tol;
 
      for (int i = 0; i < nb_ale; i++)
        {
          const int j = static_cast<int>(closest->getIJ(i, 0));
          double dist2 = 0.;
          for (int d = 0; d < dim; d++)
            {
              const double diff = ale_coords(i, d) - ext_coords(j, d);
              dist2 += diff * diff;
            }
          if (dist2 <= tol2)
            external_vertex_id_for_ale_vertex_[i] = j;
          else
            unmatched++;
        }
    }
 
  const int total_unmatched = static_cast<int>(mp_sum(unmatched));
  if (total_unmatched > 0)
    Cerr << "Domaine_ALE: " << total_unmatched
         << " ALE nodes have no external match (tol=" << tol << ")." << finl;
#endif
}
 
 
void Domaine_ALE::compute_vertex_velocities(double temps, const Probleme_base& pb)
{
  if (lagrangian_velocity_imposed_)
    {
#ifdef MEDCOUPLING_
      const Probleme_base& pb_v = ref_cast(Probleme_base, interprete().objet(external_pb_name_));
      const DoubleTab& valeurs  = pb_v.get_champ(external_velocity_field_name_).valeurs();
      ALE_mesh_velocity = 0.;
      for (int i = 0; i < nb_som(); i++)
        {
          const int j = external_vertex_id_for_ale_vertex_[i];
          if (j < 0) continue;
          for (int d = 0; d < dimension; d++)
            ALE_mesh_velocity(i, d) = valeurs(j, d);
        }
      ALE_mesh_velocity.echange_espace_virtuel();
      mesh_update_required_ = (ALE_mesh_velocity.mp_max_abs_vect() >= 1.e-12);
      return;
#else
      Cerr << "Domaine_ALE::compute_vertex_velocities: MEDCoupling required." << finl;
      Process::exit();
#endif
    }
 
  DoubleTab vit_bords(ALE_mesh_velocity);
  remplir_vitesse_bord(pb, temps, vit_bords);
  mesh_update_required_ = (vit_bords.mp_max_abs_vect() >= 1.e-12);
  if (!mesh_update_required_ && pb.schema_temps().nb_pas_dt() > 1) return;
 
  switch (meshMotionModel_)
    {
    case -1: // Boundary-only: interior nodes stay fixed
      ALE_mesh_velocity = vit_bords;
      break;
    case 0:  // Laplacian propagation
      laplacien(pb, vit_bords);
      break;
    case 1:  // Structural dynamics propagation
      structural_dynamics(pb, vit_bords, temps);
      break;
    default:
      Cerr << "Unknown ALE mesh motion model: " << meshMotionModel_ << finl;
      Process::exit();
    }
 
  Debog::verifier("Domaine_ALE::compute_vertex_velocities", ALE_mesh_velocity);
}
 
void Domaine_ALE::remplir_vitesse_bord(const Probleme_base& pb, double temps,
                                       DoubleTab& vit_bords)
{
  const int N_som = nb_som_tot();
  vit_bords = 0.;
 
  for (int n = 0; n < nb_bords_ALE; n++)
    {
      const Nom& le_nom = les_bords_ALE(n).le_nom();
      const int rang    = rang_frontiere(le_nom);
      les_champs_front[n]->associer_fr_dis_base(pb.domaine_dis().frontiere_dis(rang));
 
      DoubleTab vit_bord_ale;
      const Nom& type = les_champs_front[n]->que_suis_je();
 
      if (type == "Champ_front_ALE")
        {
          auto& c = ref_cast(Champ_front_ALE, les_champs_front[n].valeur());
          c.remplir_vit_som_bord_ALE(temps);
          vit_bord_ale = c.get_vit_som_bord_ALE();
        }
      else if (type == "Champ_front_ALE_lag")
        {
          auto& c = ref_cast(Champ_front_ALE_lag, les_champs_front[n].valeur());
          c.remplir_vit_som_bord_ALE(temps);
          vit_bord_ale = c.get_vit_som_bord_ALE();
        }
      else if (type == "Ch_front_input_ALE")
        {
          auto& c = ref_cast(Ch_front_input_ALE, les_champs_front[n].valeur());
          c.remplir_vit_som_bord_ALE(temps);
          vit_bord_ale = c.get_vit_som_bord_ALE();
        }
      else if (type == "Champ_front_ALE_Beam")
        {
          if (!eq) eq = pb.equation(0);
          auto& c = ref_cast(Champ_front_ALE_Beam, les_champs_front[n].valeur());
          c.remplir_vit_som_bord_ALE(temps);
          vit_bord_ale = c.get_vit_som_bord_ALE();
        }
      else
        {
          Cerr << "Unknown Champ_front_base type for moving boundary "
               << le_nom << ": " << type << finl;
          Process::exit();
        }
 
      // Accumulate: if multiple boundaries share nodes, only the first non-zero value is kept
      for (int dim = 0; dim < dimension; dim++)
        for (int i = 0; i < N_som; i++)
          if (les_champs_front.size() == 1 || vit_bords(i, dim) == 0.)
            vit_bords(i, dim) += vit_bord_ale(i, dim);
    }
 
  vit_bords.echange_espace_virtuel();
}
 
void Domaine_ALE::update_vertex_coordinates(double dt)
{
  // Explicit Euler update: X_{n+1} = X_n + dt * v_mesh
  sommets_ = coord_old_;
  const int N_som = nb_som_tot();
  for (int i = 0; i < N_som; i++)
    for (int d = 0; d < dimension; d++)
      coord(i, d) += ALE_mesh_velocity(i, d) * dt;
}
 
//////////////////////////////////////////////////////////////////////////////
// Laplacian mesh motion
//////////////////////////////////////////////////////////////////////////////
 
void Domaine_ALE::laplacien(const Probleme_base& pb, const DoubleTab& vit_bords)
{
  const Domaine_VF& dom_VF      = ref_cast(Domaine_VF, pb.domaine_dis());
  const Domaine_Cl_dis_base& cl = pb.equation(0).domaine_Cl_dis();
  fill_laplacian_matrix(dom_VF, cl);
  laplacien_compute_vertex_velocity(dom_VF, cl, vit_bords);
  Debog::verifier("Domaine_ALE::laplacien", ALE_mesh_velocity);
}
 
void Domaine_ALE::build_laplacian_matrix(const Domaine_VF& le_dom_VF,
                                         const Domaine_Cl_dis_base& dom_cl)
{
  const int ne_tot     = nb_elem_tot();
  const int nb_som_ele = nb_som_elem();
  const IntTab& e_s    = les_elems();
 
  // Build sparsity pattern from vertex pairs sharing an element
  Stencil stencil(0, 2);
  for (int e = 0; e < ne_tot; e++)
    for (int is = 0; is < nb_som_ele; is++)
      for (int js = is + 1; js < nb_som_ele; js++)
        {
          int i = e_s(e, is), j = e_s(e, js);
          if (i > j) { int tmp = i; i = j; j = tmp; }
          stencil.append_line(i, j);
          stencil.append_line(j, i);
          stencil.append_line(i, i);
          stencil.append_line(j, j);
        }
 
  tableau_trier_retirer_doublons(stencil);
  Matrice_Morse mat;
  Matrix_tools::allocate_morse_matrix(nb_som_tot(), nb_som_tot(), stencil, mat);
  laplacian_matrix_ = Matrice_Morse_Sym(mat);
  laplacian_matrix_.set_est_definie(1);
  Cout << "Domaine_ALE::build_laplacian_matrix: " << laplacian_matrix_.nb_lignes()
       << " rows, " << laplacian_matrix_.nb_colonnes() << " columns." << finl;
}
 
void Domaine_ALE::fill_laplacian_matrix(const Domaine_VF& le_dom_VF,
                                        const Domaine_Cl_dis_base& dom_cl)
{
  if (laplacian_matrix_.nb_lignes() != nb_som_tot())
    build_laplacian_matrix(le_dom_VF, dom_cl);
 
  laplacian_matrix_.get_set_coeff() = 0.0;
 
  const int ne_tot     = nb_elem_tot();
  const int nb_som_ele = nb_som_elem();
  const IntTab& e_f    = le_dom_VF.elem_faces();
  const IntTab& e_s    = les_elems();
  const DoubleTab& normales = le_dom_VF.face_normales();
  // In P1 vertex-based FEM on simplices, the gradient of shape function phi_i (associated
  // with node i, opposite to face i) is: grad(phi_i) = -n_i / (d * |K|)
  // where n_i is the outward face normal (with magnitude = face area) and d = space dimension.
  // The stiffness integral gives: int_K grad(phi_i).grad(phi_j) dK = (n_i.n_j) / (d^2 * |K|)
  // hence mijK = 1/d^2: 1/4 in 2D (triangles), 1/9 in 3D (tetrahedra).
  const double mijK = (dimension == 2) ? 1. / 4. : 1. / 9.;
 
  // Assemble variable-diffusivity Laplacian stiffness matrix element by element.
  // The diffusivity mu(K) = 1/|K| makes smaller cells stiffer, concentrating mesh
  // deformation in larger cells and preserving quality near walls and refined regions.
  // The element contribution is: coeffij = mu(K) * (n_i . n_j) * sign_i * sign_j * mijK / |K|
  //                                       = (n_i . n_j) * sign_i * sign_j * mijK / |K|^2
  for (int e = 0; e < ne_tot; e++)
    {
      const double volume = le_dom_VF.volumes(e);
      for (int isom = 0; isom < nb_som_ele; isom++)
        {
          int facei = e_f(e, isom), ii = e_s(e, isom);
          for (int jsom = isom + 1; jsom < nb_som_ele; jsom++)
            {
              int i = ii, facej = e_f(e, jsom), j = e_s(e, jsom);
              if (i > j) { int tmp = i; i = j; j = tmp; }
 
              double coeffij = 0.0;
              for (int d = 0; d < dimension; d++)
                coeffij += normales(facei, d) * normales(facej, d);
              coeffij *= le_dom_VF.oriente_normale(facei, e) * le_dom_VF.oriente_normale(facej, e);
              coeffij *= mijK / (volume * volume); // mu = 1/volume absorbed: mijK / |K|^2
 
              laplacian_matrix_(i, j) += coeffij;
              laplacian_matrix_(i, i) -= coeffij;
              laplacian_matrix_(j, j) -= coeffij;
            }
        }
    }
 
  // Enforce Dirichlet-like condition on all non-Neumann, non-periodic boundaries via penalization.
  // Periodic boundary nodes are left free: the Laplacian naturally interpolates the correct mesh
  // velocity for them. Penalizing them to zero would force zero mesh displacement in all directions,
  // causing mesh tearing when the domain moves (e.g. beam in x with periodic BC in z).
  for (int n_bord = 0; n_bord < le_dom_VF.nb_front_Cl(); n_bord++)
    {
      const Cond_lim& la_cl = dom_cl.les_conditions_limites(n_bord);
      const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
      if (!projection_manager_.is_neumann_boundary(le_bord.le_nom())
          && !sub_type(Periodique, la_cl.valeur()))
        {
          const int num1 = le_bord.num_premiere_face();
          const int num2 = num1 + le_bord.nb_faces();
          for (int face = num1; face < num2; face++)
            for (int isom = 0; isom < dimension; isom++)
              laplacian_matrix_(le_dom_VF.face_sommets(face, isom),
                                le_dom_VF.face_sommets(face, isom)) = ALE_PENALIZATION;
        }
    }
 
  // Also penalize nodes on internal ALE boundaries (moving solid interfaces, except Neumann)
  for (int n = 0; n < nb_bords_ALE; n++)
    {
      const Nom& le_nom = les_bords_ALE(n).le_nom();
      const int rang    = rang_frontiere(le_nom);
      if (rang < faces_bord().size() && projection_manager_.is_neumann_boundary(le_nom))
        continue;
 
      const IntTab& sommets = frontiere(le_nom).faces().les_sommets();
      const int nf   = sommets.dimension(0);
      const int nspf = sommets.dimension(1);
      for (int i = 0; i < nf; i++)
        for (int k = 0; k < nspf; k++)
          laplacian_matrix_(sommets(i, k), sommets(i, k)) = ALE_PENALIZATION;
    }
}
 
void Domaine_ALE::laplacien_compute_vertex_velocity(const Domaine_VF& le_dom_VF,
                                                    const Domaine_Cl_dis_base& dom_cl,
                                                    const DoubleTab& vit_bords)
{
  DoubleVect solution(laplacian_rhs_);
  const int nb_comp = ALE_mesh_velocity.dimension(1);
  const int nbsom   = nb_som();
 
  for (int comp = 0; comp < nb_comp; comp++)
    {
      laplacian_rhs_ = 0.;
      solution       = 0.;
 
      // Build RHS: penalized boundary values
      for (int n = 0; n < nb_bords_ALE; n++)
        {
          const Nom& le_nom = les_bords_ALE(n).le_nom();
          const int rang    = rang_frontiere(le_nom);
          if (rang < faces_bord().size() && projection_manager_.is_neumann_boundary(le_nom))
            continue;
 
          const IntTab& sommets = frontiere(le_nom).faces().les_sommets();
          const int nf   = sommets.dimension(0);
          const int nspf = sommets.dimension(1);
          for (int i = 0; i < nf; i++)
            for (int k = 0; k < nspf; k++)
              laplacian_rhs_(sommets(i, k)) = ALE_PENALIZATION * vit_bords(sommets(i, k), comp);
        }
 
      laplacian_rhs_.echange_espace_virtuel();
      Debog::verifier("Domaine_ALE::laplacien - secmem", laplacian_rhs_);
 
      if (laplacian_rhs_.mp_max_abs_vect() >= 1.e-15)
        {
          solv.resoudre_systeme(laplacian_matrix_, laplacian_rhs_, solution);
          for (int som = 0; som < nbsom; som++)
            ALE_mesh_velocity(som, comp) = solution(som);
        }
      else
        for (int som = 0; som < nbsom; som++)
          ALE_mesh_velocity(som, comp) = 0.0;
    }
 
  // Enforce exact periodicity of the mesh velocity on periodic boundary node pairs.
  // Without this, the Laplacian gives slightly different values on the two sides of each
  // periodic pair (different ghost-element stencils on each process), causing mesh-coordinate
  // drift and breaking the periodicity check in the fluid operators.
  {
    const int nb_faces_tot = le_dom_VF.face_sommets().dimension_tot(0);
    const int nb_som_face  = le_dom_VF.nb_som_face();
    IntVect fait(nb_faces_tot);
    fait = 0;
    for (int n_bord = 0; n_bord < le_dom_VF.nb_front_Cl(); n_bord++)
      {
        const Cond_lim& la_cl = dom_cl.les_conditions_limites(n_bord);
        if (!sub_type(Periodique, la_cl.valeur())) continue;
        const Periodique& la_cl_period = ref_cast(Periodique, la_cl.valeur());
        const Front_VF& le_bord = ref_cast(Front_VF, la_cl->frontiere_dis());
        const int nfin = le_bord.nb_faces_tot();
        for (int num_face = 0; num_face < nfin; num_face++)
          {
            int face = le_bord.num_face(num_face);
            if (fait(face) == 0)
              {
                int faassociee = le_bord.num_face(la_cl_period.face_associee(num_face));
                fait(face)      = 1;
                fait(faassociee) = 1;
                for (int isom = 0; isom < nb_som_face; isom++)
                  {
                    int node_src  = le_dom_VF.face_sommets(face,      isom);
                    int node_dst  = le_dom_VF.face_sommets(faassociee, isom);
                    for (int comp = 0; comp < nb_comp; comp++)
                      ALE_mesh_velocity(node_dst, comp) = ALE_mesh_velocity(node_src, comp);
                  }
              }
          }
      }
  }
 
  ALE_mesh_velocity.echange_espace_virtuel();
}
 
//////////////////////////////////////////////////////////////////////////////
// Structural dynamics mesh motion
//////////////////////////////////////////////////////////////////////////////
 
void Domaine_ALE::structural_dynamics(const Probleme_base& pb, const DoubleTab& vit_bords,
                                      double temps)
{
  const Domaine_VEF& domaine_VEF     = ref_cast(Domaine_VEF, pb.domaine_dis());
  const Domaine_Cl_VEF& domaine_Cl_VEF = ref_cast(Domaine_Cl_VEF,
                                                  pb.equation(0).domaine_Cl_dis());
  const MD_Vector& md = md_vector_sommets();
 
  // Tag nodes with imposed velocity (all nodes on any boundary face)
  IntVect tag_nodes_bords(nb_som());
  MD_Vector_tools::creer_tableau_distribue(md, tag_nodes_bords);
  tag_nodes_bords = 0;
 
  for (int n_bord = 0; n_bord < domaine_VEF.nb_front_Cl(); n_bord++)
    {
      const Front_VF& le_bord = ref_cast(Front_VF,
                                         domaine_Cl_VEF.les_conditions_limites(n_bord)->frontiere_dis());
      const int num1 = le_bord.num_premiere_face();
      const int num2 = num1 + le_bord.nb_faces();
      for (int face = num1; face < num2; face++)
        for (int isom = 0; isom < dimension; isom++)
          tag_nodes_bords[domaine_VEF.face_sommets(face, isom)] = 1;
    }
  tag_nodes_bords.echange_espace_virtuel();
 
  // Pass current coordinates to the structural model
  for (int i = 0; i < nb_som(); i++)
    for (int k = 0; k < dimension; k++)
      str_mesh_model.x(i, k) = coord(i, k);
 
  // Solve explicit dynamic problem: outputs ALE_mesh_velocity
  str_mesh_model.solveDynamicMeshProblem(
    temps, vit_bords, tag_nodes_bords, ALE_mesh_velocity,
    nb_som(), nb_elem(), nb_som_elem(), les_elems(),
    domaine_VEF.nb_faces(), domaine_VEF.nb_som_face(),
    domaine_VEF.face_sommets(), *this);
 
  // Enable grid acceleration for subsequent ICoCo sub-iterations
  if (str_mesh_model.iCoCoImplicitIteration == 0)
    str_mesh_model.acceleratedSolutionEnabled = true;
}
 
//////////////////////////////////////////////////////////////////////////////
// Post-processing fields (creer_champ / get_champ / has_champ)
//////////////////////////////////////////////////////////////////////////////
 
void Domaine_ALE::creer_champ(const Motcle& motlu, const Probleme_base& pb)
{
  Domaine::creer_champ(motlu, pb);
  const double t = pb.schema_temps().temps_courant();
 
  auto disc = [&](const Nom& type, const Nom& name, const Nom& unit,
                  int nbComp, OWN_PTR(Champ_Inc_base)& field)
  {
    pb.discretisation().discretiser_champ(type, pb.domaine_dis(), name, unit, nbComp, 1, t, field);
    champs_compris_.ajoute_champ(field);
    field->add_synonymous(name);
  };
 
  if (motlu == "ALEMeshVelocity" && !ALEMeshVelocity_)
    disc("vitesse", "ALEMeshVelocity", "m/s", dimension, ALEMeshVelocity_);
  else if (motlu == "ALEMeshTotalDisplacement" && !ALEMeshTotalDisplacement_)
    disc("vitesse", "ALEMeshTotalDisplacement", "m", dimension, ALEMeshTotalDisplacement_);
  else if (motlu == "ALEMeshStructuralPressure" && !ALEMeshStructuralPressure_)
    { assert(getMeshMotionModel() == 1); disc("champ_elem", "ALEMeshStructuralPressure", "Pa", 1, ALEMeshStructuralPressure_); }
  else if (motlu == "ALEMeshStructuralVonMises" && !ALEMeshStructuralVonMises_)
    { assert(getMeshMotionModel() == 1); disc("champ_elem", "ALEMeshStructuralVonMises", "Pa", 1, ALEMeshStructuralVonMises_); }
  else if (motlu == "ALEMeshStructuralForces" && !ALEMeshStructuralForces_)
    {
      // Note: "vitesse" discretization is used here rather than "champ_sommets"
      // to avoid an error with unknown "Champ_P1_VEF" in VEF field discretization.
      assert(getMeshMotionModel() == 1);
      disc("vitesse", "ALEMeshStructuralForces", "N", dimension, ALEMeshStructuralForces_);
    }
}
 
const Champ_base& Domaine_ALE::get_champ(const Motcle& un_nom) const
{
  if (un_nom == "ALEMeshVelocity"          && ALEMeshVelocity_)          return ALEMeshVelocity_.valeur();
  if (un_nom == "ALEMeshTotalDisplacement" && ALEMeshTotalDisplacement_) return ALEMeshTotalDisplacement_.valeur();
  if (un_nom == "ALEMeshStructuralPressure"&& ALEMeshStructuralPressure_) return ALEMeshStructuralPressure_.valeur();
  if (un_nom == "ALEMeshStructuralVonMises"&& ALEMeshStructuralVonMises_) return ALEMeshStructuralVonMises_.valeur();
  if (un_nom == "ALEMeshStructuralForces"  && ALEMeshStructuralForces_)   return ALEMeshStructuralForces_.valeur();
  throw std::runtime_error(std::string("Field ") + un_nom.getString() + " not found!");
}
 
bool Domaine_ALE::has_champ(const Motcle& un_nom, OBS_PTR(Champ_base)& ref_champ) const
{
  if (un_nom == "ALEMeshVelocity"          && ALEMeshVelocity_)          { ref_champ = ALEMeshVelocity_.valeur();          return true; }
  if (un_nom == "ALEMeshTotalDisplacement" && ALEMeshTotalDisplacement_) { ref_champ = ALEMeshTotalDisplacement_.valeur(); return true; }
  if (un_nom == "ALEMeshStructuralPressure"&& ALEMeshStructuralPressure_) { ref_champ = ALEMeshStructuralPressure_.valeur(); return true; }
  if (un_nom == "ALEMeshStructuralVonMises"&& ALEMeshStructuralVonMises_) { ref_champ = ALEMeshStructuralVonMises_.valeur(); return true; }
  if (un_nom == "ALEMeshStructuralForces"  && ALEMeshStructuralForces_)  { ref_champ = ALEMeshStructuralForces_.valeur();  return true; }
  return false;
}
 
bool Domaine_ALE::has_champ(const Motcle& un_nom) const
{
  OBS_PTR(Champ_base) dummy;
  return has_champ(un_nom, dummy);
}
 
void Domaine_ALE::get_noms_champs_postraitables(Noms& nom, Option opt) const
{
  Noms noms_compris = champs_compris_.liste_noms_compris();
  if (ALEMeshVelocity_)           noms_compris.add("ALEMeshVelocity");
  if (ALEMeshTotalDisplacement_)  noms_compris.add("ALEMeshTotalDisplacement");
  if (ALEMeshStructuralPressure_) noms_compris.add("ALEMeshStructuralPressure");
  if (ALEMeshStructuralVonMises_) noms_compris.add("ALEMeshStructuralVonMises");
  if (ALEMeshStructuralForces_)   noms_compris.add("ALEMeshStructuralForces");
 
  if (opt == DESCRIPTION)
    Cerr << que_suis_je() << " : " << noms_compris << finl;
  else
    nom.add(noms_compris);
}
 
 
void Domaine_ALE::update_coord_dom_extrait_surface()
{
  // Update vertex coordinates of any sub-domain extracted from a deformable moving boundary
  Interprete_bloc& interp = Interprete_bloc::interprete_courant();
  const Noms& noms = interp.les_noms();
  for (int i = 0; i < noms.size(); i++)
    {
      if (strcmp(interprete().objet(noms[i]).le_type(), "Domaine") != 0) continue;
      Domaine& dom_new = ref_cast(Domaine, interprete().objet(noms[i]));
      if (!sub_type(Domaine_ALE, dom_new)) continue;
      Domaine_ALE& dom_new_ale = ref_cast(Domaine_ALE, dom_new);
      if (!dom_new_ale.extrait_surf_dom_deformable_) continue;
 
      Scatter::uninit_sequential_domain(dom_new_ale);
      dom_new_ale.les_sommets() = coord_sommets();
      dom_new_ale.les_elems()   = dom_new_ale.les_elems_extrait_surf_reference();
      NettoieNoeuds::nettoie(dom_new_ale);
    }
}