Domaine_VDF.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 <Domaine_VDF.h>
#include <Domaine_Cl_VDF.h>
#include <Rectangle.h>
#include <Hexaedre.h>
#include <Faces_VDF.h>
#include <Aretes.h>
#include <Periodique.h>
#include <Dirichlet_entree_fluide_leaves.h>
#include <Neumann_sortie_libre.h>
#include <EcrFicCollecte.h>
#include <math.h>
#include <TRUSTTrav.h>
 
Implemente_instanciable(Domaine_VDF,"Domaine_VDF",Domaine_VF);
 
//#define SORT_POUR_DEBOG
 
Sortie& Domaine_VDF::printOn(Sortie& os) const
{
  Domaine_VF::printOn(os);
  os << "orientation_" << orientation_ << finl;
  os << "nb_faces_X_" << nb_faces_X_ << finl;
  os << "nb_faces_Y_" << nb_faces_Y_ << finl;
  os << "nb_faces_Z_" << nb_faces_Z_ << finl;
  os << "nb_aretes_" << nb_aretes_ << finl;
  os << "nb_aretes_joint_" << nb_aretes_joint_ << finl;
  os << "nb_aretes_coin_" << nb_aretes_coin_ << finl;
  os << "nb_aretes_bord_" << nb_aretes_bord_ << finl;
  os << "nb_aretes_mixtes_" << nb_aretes_mixtes_ << finl;
  os << "nb_aretes_internes_" << nb_aretes_internes_ << finl;
  Qdm_.ecrit(os << "Qdm : ");
  os << h_x_ << " " << h_y_ << " " << h_z_ << finl;
  return os;
}
 
Entree& Domaine_VDF::readOn(Entree& is)
{
  Domaine_VF::readOn(is);
  is >> orientation_;
  is >> nb_faces_X_;
  is >> nb_faces_Y_;
  is >> nb_faces_Z_;
  is >> nb_aretes_;
  is >> nb_aretes_joint_;
  is >> nb_aretes_coin_;
  is >> nb_aretes_bord_;
  is >> nb_aretes_mixtes_;
  is >> nb_aretes_internes_;
  Qdm_.lit(is);
  is >> h_x_ >> h_y_ >> h_z_;
  return is;
}
/*! @brief Returns a newly allocated Faces_VDF object.
 *
 * @return Pointer to the newly created Faces_VDF instance.
 */
Faces* Domaine_VDF::creer_faces()
{
  Faces_VDF* les_faces_vdf=new(Faces_VDF);
  return les_faces_vdf;
}
 
/*! Override to do nothing. Not necessary for VDF.
 */
void Domaine_VDF::prepare_elem_non_std(Faces&)
{
}
 
/*! @brief Override. Compute sorting key so that internal faces are sorted by their orientation first (X, Y, Z)
 */
void Domaine_VDF::compute_sort_key(Faces& les_faces, IntTab& sort_key)
{
  // Compute the orientation of real faces
  Faces_VDF& les_faces_vdf=ref_cast(Faces_VDF, les_faces);
  les_faces_vdf.calculer_orientation(orientation_, nb_faces_X_, nb_faces_Y_, nb_faces_Z_);
 
  const int nb_faces_front = domaine().nb_faces_frontiere();
 
  // Build an integer key for sorting the faces:
  //  orientation * nb_faces + face_index
  // Sorting by ascending order of this key sorts by orientation
  // while preserving the original order of faces with the same orientation
  const int nb_faces = les_faces_vdf.nb_faces();
 
  nb_faces_std_ = 0;
  sort_key.resize(nb_faces, 2);
  // Do not sort the boundary faces, which remain at the beginning:
  for (int i = 0; i < nb_faces_front; i++)
    {
      sort_key(i, 0) = i;
      sort_key(i, 1) = i;
    }
 
  for (int i=nb_faces_front; i < nb_faces; i++)
    {
      const int ori = orientation_[i];
      sort_key(i, 0) = ori * nb_faces + i;
      sort_key(i, 1) = i;
    }
  nb_faces_std_ = nb_faces - nb_faces_front;  // not as funny as in VEF ... :-)
}
 
/*! @brief Override to also renumber orientation_ member.
 */
void Domaine_VDF::renumber_faces(Faces& les_faces, IntTab& sort_key)
{
  Domaine_VF::renumber_faces(les_faces, sort_key);
 
  const int nb_faces_front = domaine().nb_faces_frontiere();
  const int nb_faces = les_faces.nb_faces();
 
  IntVect old_orien(orientation_);
  for (int i=nb_faces_front; i < nb_faces; i++)
    {
      const int idx = sort_key(i, 1);
      orientation_[i] = old_orien[idx];
    }
}
 
 
/*! @brief Calls Domaine_VF::discretiser(), computes element gravity centers,
 *
 * @brief fills element/face connectivities, reorders face/element connectivities
 *  so that el1 is on the left and el2 on the right,
 *  computes volumetric and surface porosities,
 *  generates edges, and computes mesh steps.
 *
 */
void Domaine_VDF::discretiser()
{
  Domaine_VF::discretiser();
 
  Domaine& domaine_geom=domaine();
 
  // Check the consistency between the geometric element and the discretization
 
  const Elem_geom_base& elem_geom = domaine_geom.type_elem().valeur();
 
  if (dimension == 2)
    {
      if (!sub_type(Rectangle,elem_geom))
        {
          Cerr << " The geometric element type " << elem_geom.que_suis_je() << " is incorrect" << finl;
          Cerr << " Only the Rectangle and Rectangle_Axi types are compatible with VDF discretization in dimension 2" << finl;
          exit();
        }
    }
  else if (dimension == 3)
    {
      if (!sub_type(Hexaedre,elem_geom))
        {
          Cerr << " The geometric element type " << elem_geom.que_suis_je() << " is incorrect" << finl;
          Cerr << " Only the Hexaedre and Hexaedre_Axi types are compatible with VDF discretization in dimension 3" << finl;
          exit();
        }
    }
  // Compute the orientation of virtual faces
  // Note BM: for compatibility with version 1.5, orientation_
  //  has no virtual space.
  MD_Vector md_nul;
  creer_tableau_faces(orientation_);
  orientation_.echange_espace_virtuel();
  orientation_.set_md_vector(md_nul); // Detach the parallel structure
 
  // Apply the VDF convention on face_voisin:
  // The element face_voisin(i,0) must have an "ori" coordinate smaller than the face
  // and the element face_voisin(i,1) must have a larger coordinate.
  {
    const int nbr_faces_tot = face_voisins_.dimension_tot(0);
    for (int i_face = 0; i_face < nbr_faces_tot; i_face++)
      {
        const int elem0 = face_voisins_(i_face, 0);
        const int elem1 = face_voisins_(i_face, 1);
        const int ori = orientation_[i_face];
        const double x_face = xv(i_face, ori);
        double delta = 0.;
        // Element 0 must have an "ori" coordinate smaller than the face
        // and element 1 must have a larger coordinate.
        if (elem0 >= 0)
          {
            delta = x_face - xp(elem0, ori);
          }
        else
          {
            assert(elem1 >= 0); // Otherwise, critical error: neighboring elements not filled in
            delta = xp(elem1, ori) - x_face;
          }
        if (delta < 0)
          {
            // Swap the two elements
            face_voisins_(i_face, 0) = elem1;
            face_voisins_(i_face, 1) = elem0;
          }
      }
  }
 
  calculer_volumes_entrelaces();
  genere_aretes();
  calcul_h();
  remplir_face_normales();
  Cerr << "The Domaine_VDF object has been filled successfully " << finl;
 
}
 
void Domaine_VDF::remplir_face_normales()
{
  // Fill the face_normales_ array;
  //  Note: the face_voisins array is not exactly a distributed array. A face does not have its two neighbors in the
  //  same order on all processors that own the face.
  //  Therefore the face normal can change direction from one processor to another, including for joint faces.
  face_normales_ = xv_; // already has // structure
  face_normales_ = 0.;
 
  for (int f = 0; f < nb_faces_tot(); f++)
    {
      int ori = orientation(f);
      face_normales_(f,ori) = face_surfaces(f);
    }
}
 
/*! @brief Fills the staggered volumes (volumes entrelaces).
 *
 */
void Domaine_VDF::calculer_volumes_entrelaces()
{
  Cerr << "Computing staggered volumes" << finl;
  creer_tableau_faces(volumes_entrelaces_);
  volumes_entrelaces_dir_.resize(nb_faces(), 2);
  creer_tableau_faces(volumes_entrelaces_dir_);
 
  const int nb_faces_front = premiere_face_int();
  const int nbf = nb_faces();
  for (int num_face = 0; num_face<nbf; num_face++)
    {
      const double f = (num_face < nb_faces_front) ? 2. : 1.;
      for (int dir = 0; dir < 2; dir ++)
        {
          int elem = face_voisins_(num_face, dir);
          if (elem != -1)
            {
              if ((axi) && (orientation_[num_face]==0))
                volumes_entrelaces_dir_(num_face, dir) = f * 0.5 * xv_(num_face, 0) * volumes(elem) / xp(elem, 0);
              else if ((bidim_axi) && (orientation_[num_face]==0))
                volumes_entrelaces_dir_(num_face, dir) = volume_entrelace_axi(std::fabs(xv(num_face, 0)), std::fabs(xp(elem, 0)), dim_elem(elem, 1));
              else
                volumes_entrelaces_dir_(num_face, dir) = f * 0.5 * volumes(elem);
 
              volumes_entrelaces_[num_face] += volumes_entrelaces_dir_(num_face, dir);
            }
        }
    }
  volumes_entrelaces_.echange_espace_virtuel();
  volumes_entrelaces_dir_.echange_espace_virtuel();
}
 
static inline int face_vois(const Domaine_VDF& zvdf, const Domaine& domaine, int face, int i)
{
  int elem=zvdf.face_voisins(face,i);
  if(elem==-1)
    {
      face=domaine.face_bords_interne_conjuguee(face);
      if(face!=-1)
        {
          elem=zvdf.face_voisins(face,i);
          assert(elem!=-1);
        }
    }
  return elem;
}
 
/*! @brief Generates the edge connectivity arrays.
 *
 */
void Domaine_VDF::genere_aretes()
{
  Cerr << "Generating edges" << finl;
 
  Domaine& mon_dom=domaine();
  //  int nb_poly = le_dom.nb_elem();
  int nb_poly_tot = mon_dom.nb_elem_tot();
  // Estimate with upper bound on the number of edges: nb_aretes_plus
  nb_aretes_=-1; // cf Aretes::affecter
  int nb_aretes_plus=-1;
  if (dimension == 2)
    nb_aretes_plus = mon_dom.nb_som();
  else if (dimension == 3)
    nb_aretes_plus = mon_dom.nb_som()*3;
  Cerr << "Creating edges " << finl;
  Aretes les_aretes(nb_aretes_plus);
 
  // Loop over all elements:
  int el1, el2, el3, el4;
  int face12, face13, face34, face24;
 
  int nb_dir;
  if(dimension==2) nb_dir=1;
  else nb_dir=3;
  IntVect gauche(nb_dir);
  gauche(0)=0;
  if(dimension==3)
    {
      gauche(1)=0;
      gauche(2)=1;
    }
  IntVect droite(nb_dir);
  droite(0)=dimension;
  if(dimension==3)
    {
      droite(1)=dimension;
      droite(2)=dimension+1;
    }
  IntVect haut(nb_dir);
  haut(0)=dimension+1;
  if(dimension==3)
    {
      haut(1)=dimension+2;
      haut(2)=dimension+2;
    }
  IntVect bas(nb_dir);
  bas(0)=1;
  if(dimension==3)
    {
      bas(1)=2;
      bas(2)=2;
    }
 
  int coin=-1;
  int bord=0;
  int mixte=1;
  int interne=2;
  // Detection of plates (2 boundary faces overlap):
  IntVect est_une_plaque(nb_faces());
  creer_tableau_faces(est_une_plaque);
  est_une_plaque=0;
  // Loop over boundary faces
  ArrOfDouble P1(3), P2(3);
  P1=0;
  P2=0;
  double eps = 10.0*Objet_U::precision_geom;
  for (int face=0; face<premiere_face_int(); face++)
    {
      int ori = orientation(face);
      for (int i = 0; i < dimension; i++)
        {
          P1[i] = xv(face, i) + (ori == i ? eps : 0);
          P2[i] = xv(face, i) - (ori == i ? eps : 0);
        }
      int elem1 = domaine().chercher_elements(P1[0], P1[1], P1[2]);
      int elem2 = domaine().chercher_elements(P2[0], P2[1], P2[2]);
      if (elem1 >= 0 && elem2 >= 0 && elem1 != elem2)
        {
          // Find 2 points P1 and P2 in cells so:
          est_une_plaque(face) = 1;
          //Journal() << "We detect an internal boundary on face " << face << " between elements " << elem1 << " and " << elem2 << finl;
        }
    }
  est_une_plaque.echange_espace_virtuel();
 
  for(int dir=0; dir<nb_dir; dir++)
    for (el1=0; el1<nb_poly_tot; el1++)
      {
        // Generate the edge at the top-right of el1
        face12=elem_faces(el1,droite(dir));
        face13=elem_faces(el1,haut(dir));
 
        el2=face_vois(*this, mon_dom, face12, 1);
 
        if(el2>-1)
          face24=elem_faces(el2,haut(dir));
        else
          face24=-1;
 
        el3=face_vois(*this, mon_dom, face13, 1);
 
        if(el3>-1)
          face34=elem_faces(el3,droite(dir));
        else
          face34=-1;
 
        if( (el2>-1) && (el3>-1))
          el4=face_vois(*this, mon_dom, face24, 1);
        else if((el2>-1))
          el4=face_vois(*this, mon_dom, face24, 1);
        else if((el3>-1))
          el4=face_vois(*this, mon_dom, face34, 1);
        else
          el4=-1;
 
        if ( (el2==-1) && (el4>=0) )
          face24=elem_faces(el4,bas(dir));
        if ( (el3==-1) && (el4>=0) )
          face34=elem_faces(el4,gauche(dir));
 
        const int nb_f = nb_faces();
        if (el2 > -1 && el3 > -1 && el4 > -1) // internal edge
          les_aretes.affecter(nb_aretes_, dir, interne, nb_f, face13, face24, face12, face34, est_une_plaque);
        else if ( (el3 > -1 && el4 > -1) ||
                  (el2 > -1 && el4 > -1) ||
                  (el2 > -1 && el3 > -1) ) // mixed edge
          les_aretes.affecter(nb_aretes_, dir, mixte, nb_f, face13, face24, face12, face34, est_une_plaque);
        else if (el2 > -1) // boundary edge
          les_aretes.affecter(nb_aretes_, dir, bord, nb_f, face13, face24, face12, 1, est_une_plaque);
        else if (el3 > -1) // boundary edge
          les_aretes.affecter(nb_aretes_, dir, bord, nb_f, face12, face34, face13, 1, est_une_plaque);
        else // corner edge
          les_aretes.affecter(nb_aretes_, dir, coin, nb_f, face13, face24, face12, face34, est_une_plaque);
        //Cerr << "elements_haut_droit " << el1 << " " << el2 << " " << el3 << " " << el4 << finl;
        //Journal() << "Provisoire faces arete: " << face12 << " " << face13 << " " << face24 << " " << face34 << finl;
 
        // For corners or boundaries:
        face13 = elem_faces(el1, bas(dir));
        el3 = face_vois(*this, mon_dom, face13, 0);
        if (el3 < 0) // Generate the edge at the bottom-right of el1
          // if the cell below el1 does not exist
          {
            if (el2 >= 0)
              {
                face24 = elem_faces(el2, bas(dir));
                el4 = face_vois(*this, mon_dom, face24, 0);
                if (el4 < 0) // boundary edge
                  les_aretes.affecter(nb_aretes_, dir, bord, nb_f, face13, face24, face12, -1, est_une_plaque);
                else // mixed edge
                  {
                    face34 = elem_faces(el4, gauche(dir));
                    if (face_vois(*this, mon_dom, face34, 0) == -1)
                      les_aretes.affecter(nb_aretes_, dir, mixte, nb_f, face13, face24, face34, face12, est_une_plaque);
                  }
              }
            else   // corner edge
              les_aretes.affecter(nb_aretes_, dir, coin, nb_f, face13, -1, -1, face12, est_une_plaque);
            //Cerr << "elements_bas_droit " << el1 << " " << el2 << " " << el3 << " " << el4 << finl;
          }
 
 
        face13 = elem_faces(el1, gauche(dir));
        el3 = face_vois(*this, mon_dom, face13, 0);
        if (el3 < 0) // Generate the edge at the top-left of el1
          // if the cell to the top-left of el1 does not exist
          {
            face12 = elem_faces(el1, haut(dir));
            el2 = face_vois(*this, mon_dom, face12, 1);
            if (el2 >= 0)
              {
                face24 = elem_faces(el2, gauche(dir));
                el4 = face_vois(*this, mon_dom, face24, 0);
                if (el4 < 0) // boundary edge
                  les_aretes.affecter(nb_aretes_, dir, bord, nb_f, face13, face24, face12, -1, est_une_plaque);
              }
            else   // corner edge
              les_aretes.affecter(nb_aretes_, dir, coin, nb_f, -1, face12, -1, face13, est_une_plaque);
            //Cerr << "elements_haut_gauche " << el1 << " " << el2 << " " << el3 << " " << el4 << finl;
          }
 
        // Generate the edge at the bottom-left of el1
        face12 = elem_faces(el1, gauche(dir));
        el2 = face_vois(*this, mon_dom, face12, 0);
        face13 = elem_faces(el1, bas(dir));
        el3 = face_vois(*this, mon_dom, face13, 0);
        if ((el2 < 0) && (el3 < 0)) // corner edge
          {
            les_aretes.affecter(nb_aretes_, dir, coin, nb_f, face13, -1, face12, -1, est_une_plaque);
            //Cerr << "elements_bas_gauche " << el1 << " " << el2 << " " << el3 << " " << el4 << finl;
          }
      }
  nb_aretes_++;
  les_aretes.dimensionner(nb_aretes_);
  // Cerr << "Tri des aretes " << finl;
  les_aretes.trier(nb_aretes_coin_,
                   nb_aretes_bord_,
                   nb_aretes_mixtes_,
                   nb_aretes_internes_);
#ifdef SORT_POUR_DEBOG
 
  les_aretes.trier_pour_debog(nb_aretes_coin_,
                              nb_aretes_bord_,
                              nb_aretes_mixtes_,
                              nb_aretes_internes_,xv());
#endif
  nb_aretes_joint_=0;
  assert(nb_aretes_==nb_aretes_coin_+nb_aretes_bord_+nb_aretes_mixtes_
         +nb_aretes_internes_+nb_aretes_joint_);
  Qdm_.ref(les_aretes.faces());
  /*
  // Loop to verify where the internal walls are
  for (int i=0; i<Qdm_.dimension(0); i++)
    {
      int nb_plaques = 0;
      for (int j = 0; j < Qdm_.dimension(1); j++)
        nb_plaques += Qdm_(i, j) >= 0 ? est_une_plaque(Qdm_(i, j)) : 0;
      if (nb_plaques > 0)
      Journal() << "Provisoire arete " << i << " a " << nb_plaques << " paroi internes." << finl;
    }
  Journal() << "Aretes coin   = " << nb_aretes_coin_ << finl;
    Journal() << "Aretes bord   = " << nb_aretes_bord_ << finl;
    Journal() << "Aretes mixte  = " << nb_aretes_mixtes_ << finl;
    Journal() << "Aretes interne= " << nb_aretes_internes_ << finl;
    Journal() << "Aretes joint  = " << nb_aretes_joint_ << finl;
     */
  //Cerr << "Qdm : " << Qdm_ << finl;
}
 
/*! @brief Computes the mesh step sizes h_x_, h_y_, h_z_.
 *
 */
void Domaine_VDF::calcul_h()
{
  Domaine& domaine_geom=domaine();
  IntVect numfa(domaine_geom.nb_faces_elem());
  const double deux_pi = M_PI * 2.0;
  const int D = dimension;
 
  for (int e = 0; e < domaine_geom.nb_elem(); e++)
    {
      for (int j = 0; j < domaine_geom.nb_faces_elem(); j++)
        numfa[j] = elem_faces(e, j);
 
      h_x_ = std::min(h_x_, xv_(numfa[D], 0) - xv_(numfa[0], 0));
      double hy_loc;
      if (axi)
        {
          double d_teta = xv_(numfa[D + 1], 1) - xv_(numfa[1], 1);
          if (d_teta < 0) d_teta += deux_pi;
          hy_loc = d_teta * xv_(numfa[1], 0);
        }
      else hy_loc = xv_(numfa[D + 1], 1) - xv_(numfa[1], 1);
      h_y_ = std::min(h_y_, hy_loc);
      if (dimension == 3) h_z_ = std::min(h_z_, xv_(numfa[5], 2) - xv_(numfa[2], 2));
    }
  h_x_ = mp_min(h_x_);
  h_y_ = mp_min(h_y_);
  h_z_ = mp_min(h_z_);
}
 
void Domaine_VDF::modifier_pour_Cl(const Conds_lim& conds_lim)
{
  // Cerr << "Domaine_VDF::Modifier_pour_Cl" << finl;
  Domaine_VF::modifier_pour_Cl(conds_lim);
 
  int fac3;
  int nb_cond_lim=conds_lim.size();
  IntTab aretes_coin_traitees(nb_aretes_coin());
  aretes_coin_traitees = -1;
 
  for (int num_cond_lim=0; num_cond_lim<nb_cond_lim; num_cond_lim++)
    {
      // for cl ...
      //Journal() << "On traite la cl num : " << num_cond_lim << finl;
      const Cond_lim_base& cl = conds_lim[num_cond_lim].valeur();
      if (sub_type(Periodique, cl))
        {
          // if cl perio
 
          // Modification of the Qdm_ array for periodic edge types
 
          const Domaine_Cl_VDF& domaine_Cl_VDF = ref_cast(Domaine_Cl_VDF,cl.domaine_Cl_dis());
 
          int ndeb_arete = premiere_arete_bord();
          int fac1,fac2,sign,num_face,elem1,n_type;
          for (int n_arete=ndeb_arete; n_arete<ndeb_arete+nb_aretes_bord(); n_arete++)
            {
              // for n_arete
              n_type=domaine_Cl_VDF.type_arete_bord(n_arete-ndeb_arete);
 
              if (n_type == TypeAreteBordVDF::PERIO_PERIO) // periodic edge type
                {
                  //if arete perio
                  fac1 = Qdm_(n_arete,0);
                  fac2 = Qdm_(n_arete,1);
                  fac3 = Qdm_(n_arete,2);
                  sign = Qdm_(n_arete,3);
 
                  const Front_VF& la_frontiere_dis = ref_cast(Front_VF,cl.frontiere_dis());
                  int ndeb = la_frontiere_dis.num_premiere_face();
                  int nfin = ndeb + la_frontiere_dis.nb_faces();
                  if ( ( ( ndeb <= fac1) && (fac1 < nfin) ) ||
                       ( ( ndeb <= fac2) && (fac2 < nfin) )
                     )
                    {
                      if (sign == 1)
                        elem1 = face_voisins(fac1,1);
                      else
                        elem1 = face_voisins(fac1,0);
                      num_face = dimension+orientation_[fac3];
 
                      if (sign == -1)
                        {
                          Qdm_(n_arete,3) = fac3;
                          Qdm_(n_arete,2) = elem_faces(elem1,num_face);
                        }
                      else
                        Qdm_(n_arete,3) = elem_faces(elem1,num_face);
                    }
                }
            }
          // handle corner edges
          int ndeb_arete_coin = premiere_arete_coin();
          Cerr << "Modifying Qdm for corner edges  num_cond_lim=" << num_cond_lim <<  finl;
          int fac4;
 
          for (int n_arete=ndeb_arete_coin; n_arete<ndeb_arete_coin+nb_aretes_coin(); n_arete++)
            {
              if (aretes_coin_traitees[n_arete] != 1)
                {
                  fac1 = Qdm_(n_arete,0);
                  fac2 = Qdm_(n_arete,1);
                  fac3 = Qdm_(n_arete,2);
                  fac4 = Qdm_(n_arete,3);
 
                  // Retrieve face indices that are not equal to -1
                  IntVect f(2);
                  f = -2;
                  int i=0;
                  if (fac1 != -1)
                    {
                      f(i) = fac1;
                      i++;
                    }
                  if (fac2 != -1)
                    {
                      f(i) = fac2;
                      i++;
                    }
                  if (fac3 != -1)
                    {
                      f(i) = fac3;
                      i++;
                    }
                  if (fac4 != -1)
                    {
                      f(i) = fac4;
                      i++;
                    }
 
                  // Check whether the face is a periodic face
 
                  const Front_VF& la_frontiere_dis = ref_cast(Front_VF,cl.frontiere_dis());
                  int ndeb = la_frontiere_dis.num_premiere_face();
                  int nfin = ndeb + la_frontiere_dis.nb_faces();
                  int  elem, fac,dim0,dim1,indic_f0=-100,indic_f1=-100;
 
                  n_type=domaine_Cl_VDF.type_arete_coin(n_arete-ndeb_arete_coin);
                  dim1 = orientation_[f(1)];
                  dim0 = orientation_[f(0)];
 
                  for (int j=0; j<2; j++)
                    for (int k=0; k<2; k++)
                      if ((face_voisins(f(0),j) == face_voisins(f(1),k))  && (face_voisins(f(0),j)!=-1) )
                        {
                          indic_f0 = j;
                          indic_f1 = k;
                        }
 
                  if ((n_type == TypeAreteCoinVDF::PERIO_PAROI) || (n_type == TypeAreteCoinVDF::PERIO_FLUIDE))// perio-wall corner edge
                    {
                      if ((f(0) >= ndeb)&&(f(0) < nfin))
                        {
                          Qdm_(n_arete,2)=f(0);
                          Qdm_(n_arete,indic_f0)=f(1);
 
                          elem = face_voisins(f(0),1-indic_f0);
                          fac = elem_faces(elem,dim1+(1-indic_f1)*dimension);
                          Qdm_(n_arete,1-indic_f0)=fac;
 
                          Qdm_(n_arete,3)=1-2*indic_f1;
 
                          //                          Cerr << "n_arete=" << n_arete << "  OK!!" << finl;
                          aretes_coin_traitees[n_arete] = 1;
                        }
                      else
                        {
                          if ((f(1) >= ndeb)&&(f(1) < nfin))
                            {
                              Qdm_(n_arete,2)=f(1);
                              Qdm_(n_arete,indic_f1)=f(0);
 
                              elem = face_voisins(f(1),1-indic_f1);
                              fac = elem_faces(elem,dim0+(1-indic_f0)*dimension);
                              Qdm_(n_arete,1-indic_f1)=fac;
 
                              Qdm_(n_arete,3)=1-2*indic_f0;
 
                              //                              Cerr << "n_arete=" << n_arete << "  OK!!" << finl;
                              aretes_coin_traitees[n_arete] = 1;
                            }
                          else
                            {
                              // Cerr<<"Attention cas non traite lors du remplissage du tableau Qdm"<<finl;
                              // exit();
                              ;
                            }
                        }
                    }
                  else if (n_type == TypeAreteCoinVDF::PERIO_PERIO) // perio-perio corner edge
                    {
                      if ((f(0) >= ndeb)&&(f(0) < nfin))
                        {
                          Qdm_(n_arete,2+indic_f1)=f(0);
                          Qdm_(n_arete,indic_f0)=f(1);
 
                          elem = face_voisins(f(0),1-indic_f0);
                          fac = elem_faces(elem,dim1+(1-indic_f1)*dimension);
                          Qdm_(n_arete,1-indic_f0)=fac;
 
                          elem = face_voisins(f(1),1-indic_f1);
                          fac = elem_faces(elem,dim0+(1-indic_f0)*dimension);
                          Qdm_(n_arete,3-indic_f1)=fac;
 
                          //                            Cerr << "n_arete=" << n_arete << "  OK!!" << finl;
                          aretes_coin_traitees[n_arete] = 1;
                        }
                      else
                        {
                          if ((f(1) >= ndeb)&&(f(1) < nfin))
                            {
                              Qdm_(n_arete,2+indic_f0)=f(1);
                              Qdm_(n_arete,indic_f1)=f(0);
 
                              elem = face_voisins(f(1),1-indic_f1);
                              fac = elem_faces(elem,dim0+(1-indic_f0)*dimension);
                              Qdm_(n_arete,1-indic_f1)=fac;
 
                              elem = face_voisins(f(0),1-indic_f0);
                              fac = elem_faces(elem,dim1+(1-indic_f1)*dimension);
                              Qdm_(n_arete,3-indic_f0)=fac;
 
                              //                                Cerr << "n_arete=" << n_arete << "  OK!!" << finl;
                              aretes_coin_traitees[n_arete] = 1;
                            }
                        }
                    }
 
                  else
                    {
                      Cerr << "Warning: cases involving other corner edge types "  << finl;
                      Cerr << "than perio-perio or perio-wall are not handled!!" << finl;
                    }
                }
            }
        }
 
      // Modif OC 01/2005 to handle wall/fluid corner types
      else if (sub_type(Dirichlet_entree_fluide, cl) || sub_type(Neumann_sortie_libre, cl) )
        {
          // if cl of wall type
 
          const Domaine_Cl_VDF& domaine_Cl_VDF = ref_cast(Domaine_Cl_VDF,cl.domaine_Cl_dis());
 
          int fac1,fac2,n_type;
          int ndeb_arete_coin = premiere_arete_coin();
          Cerr << "Modifying Qdm for corner edges touching a wall num_cond_lim=" << num_cond_lim <<  finl;
          int fac4;
 
          for (int n_arete=ndeb_arete_coin; n_arete<ndeb_arete_coin+nb_aretes_coin(); n_arete++)
            {
              if (aretes_coin_traitees[n_arete] != 1)
                {
                  fac1 = Qdm_(n_arete,0);
                  fac2 = Qdm_(n_arete,1);
                  fac3 = Qdm_(n_arete,2);
                  fac4 = Qdm_(n_arete,3);
 
                  // Retrieve face indices that are not equal to -1
                  IntVect f(2);
                  f = -2;
                  int i=0;
                  if (fac1 != -1)
                    {
                      f(i) = fac1;
                      i++;
                    }
                  if (fac2 != -1)
                    {
                      f(i) = fac2;
                      i++;
                    }
                  if (fac3 != -1)
                    {
                      f(i) = fac3;
                      i++;
                    }
                  if (fac4 != -1)
                    {
                      f(i) = fac4;
                      i++;
                    }
 
                  // Check whether the face is a wall-type face
 
                  const Front_VF& la_frontiere_dis = ref_cast(Front_VF,cl.frontiere_dis());
                  int ndeb = la_frontiere_dis.num_premiere_face();
                  int nfin = ndeb + la_frontiere_dis.nb_faces();
                  int indic_f0=-100,indic_f1=-100;
 
                  n_type=domaine_Cl_VDF.type_arete_coin(n_arete-ndeb_arete_coin);
 
                  for (int j=0; j<2; j++)
                    for (int k=0; k<2; k++)
                      if ((face_voisins(f(0),j) == face_voisins(f(1),k)) && (face_voisins(f(0),j)!=-1) )
                        {
                          indic_f0 = j;
                          indic_f1 = k;
                        }
 
                  if ((n_type == TypeAreteCoinVDF::PAROI_FLUIDE) || (n_type == TypeAreteCoinVDF::FLUIDE_PAROI) || (n_type == TypeAreteCoinVDF::FLUIDE_FLUIDE))// wall-fluid corner edge
                    {
                      if ((f(0) >= ndeb)&&(f(0) < nfin))
                        {
                          Qdm_(n_arete,0)=f(1);
                          Qdm_(n_arete,1)=f(1);
                          Qdm_(n_arete,2)=f(0);
                          Qdm_(n_arete,3)=1-2*indic_f1;
                          aretes_coin_traitees[n_arete] = 1;
                        }
                      else
                        {
                          if ((f(1) >= ndeb)&&(f(1) < nfin))
                            {
                              Qdm_(n_arete,0)=f(0);
                              Qdm_(n_arete,1)=f(0);
                              Qdm_(n_arete,2)=f(1);
                              Qdm_(n_arete,3)=1-2*indic_f0;
                              aretes_coin_traitees[n_arete] = 1;
                            }
                          else
                            {
                              // Cerr<<"Attention cas non traite lors du remplissage du tableau Qdm"<<finl;
                              // exit();
                              ;
                            }
                        }
                    }
                }
            }
        }
 
      //   Cerr << "Qdm : " << Qdm_ << finl;
      //   Cerr << "aretes_coin_traitees : " << aretes_coin_traitees << finl;
    }
  //Journal() << "domain after modification for Cl : " << *this << finl;
 
 
  // PQ : 10/10/05 : periodic faces have a double contribution;
  //                      the call to marquer_faces_double_contrib is done in this method
  //                      in order to have access to conds_lim.
 
 
  Domaine_VF::marquer_faces_double_contrib(conds_lim);
}
 
/*! @brief Fills the face_voisins_fictifs_ array. Does NOT create fictitious elements.
 *
 * @param zcldisbase The discretized boundary condition domain.
 */
 
void Domaine_VDF::creer_elements_fictifs(const Domaine_Cl_dis_base& zcldisbase)
{
  Cerr << "Domaine_VDF::creer_elements_fictifs() called" << finl;
  const Domaine_Cl_VDF& zclvdf = ref_cast(Domaine_Cl_VDF,zcldisbase);
  if (face_voisins_fictifs_.size() == 0)
    {
      face_voisins_fictifs_.resize(nb_faces_bord(),2);
      face_voisins_fictifs_ = -1;
      int compteur = nb_elem_tot();
      int face,ndeb,nfin;
 
      for (int n_bord=0; n_bord<nb_front_Cl(); n_bord++)
        {
          const Cond_lim& la_cl = zclvdf.les_conditions_limites(n_bord);
          const Front_VF& le_bord = ref_cast(Front_VF,la_cl->frontiere_dis());
          if ( (sub_type(Dirichlet_entree_fluide,la_cl.valeur())) ||
               (sub_type(Neumann_sortie_libre,la_cl.valeur())) )
            {
              ndeb = le_bord.num_premiere_face();
              nfin = ndeb + le_bord.nb_faces();
              for (face=ndeb; face<nfin; face++)
                if (face_voisins(face,0) != -1)
                  face_voisins_fictifs_(face,1) = compteur++;
                else
                  face_voisins_fictifs_(face,0) = compteur++ ;
            }
        }
    }
  Cerr << "size "<<face_voisins_fictifs_.size()-nb_elem()<<finl;
}
 
// Fills the dist array with the normal distances to the boundary
// for the boundary faces of name nom_bord.
// The dist array is resized by this method.
DoubleVect& Domaine_VDF::dist_norm_bord(DoubleVect& dist, const Nom& nom_bord) const
{
  if (axi)
    {
      for (int n_bord=0; n_bord<les_bords_.size(); n_bord++)
        {
          const Front_VF& fr_vf = front_VF(n_bord);
          if (fr_vf.le_nom() == nom_bord)
            {
              dist.resize(fr_vf.nb_faces());
              int ndeb = fr_vf.num_premiere_face();
              for (int face=ndeb; face<ndeb+fr_vf.nb_faces(); face++)
                dist(face-ndeb) = dist_norm_bord_axi(face);
            }
        }
    }
  else
    {
      for (int n_bord=0; n_bord<les_bords_.size(); n_bord++)
        {
          const Front_VF& fr_vf = front_VF(n_bord);
          if (fr_vf.le_nom() == nom_bord)
            {
              dist.resize(fr_vf.nb_faces());
              int ndeb = fr_vf.num_premiere_face();
              for (int face=ndeb; face<ndeb+fr_vf.nb_faces(); face++)
                dist(face-ndeb) = dist_norm_bord(face);
            }
        }
    }
  return dist;
}
 
void Domaine_VDF::init_virt_e_map() const
{
  IntTrav p_e(0, 2);
  domaine().creer_tableau_elements(p_e);
  for (int e = 0; e < nb_elem() ; e++) p_e(e, 0) = Process::me(), p_e(e, 1) = e;
  p_e.echange_espace_virtuel();
  for (int e = nb_elem() ; e < nb_elem_tot() ; e++) virt_e_map[ {{ p_e(e, 0), p_e(e, 1) }}] = e;
}