Domaine_VEF.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 <Connectivite_som_elem.h>
#include <Check_espace_virtuel.h>
#include <MD_Vector_composite.h>
#include <VEF_discretisation.h>
#include <EcrFicPartageBin.h>
#include <Static_Int_Lists.h>
#include <MD_Vector_tools.h>
#include <Domaine_Cl_VEF.h>
#include <Quadrangle_VEF.h>
#include <Octree_Double.h>
#include <Hexaedre_VEF.h>
#include <Domaine_VEF.h>
#include <EFichierBin.h>
#include <Quadri_VEF.h>
#include <Periodique.h>
#include <TRUSTLists.h>
#include <Periodique.h>
#include <Rectangle.h>
#include <Tetra_VEF.h>
#include <Conds_lim.h>
#include <Tetraedre.h>
#include <ArrOfBit.h>
#include <Hexaedre.h>
#include <Hexa_VEF.h>
#include <Triangle.h>
#include <Tri_VEF.h>
#include <Domaine.h>
#include <Debog.h>
 
Implemente_instanciable(Domaine_VEF, "Domaine_VEF", Domaine_VF);
 
Sortie& Domaine_VEF::ecrit(Sortie& os) const
{
  Domaine_VF::printOn(os);
  os << "____ h_carre " << finl;
  os << h_carre << finl;
  os << "____ type_elem_ " << finl;
  os << type_elem_.valeur() << finl;
  os << "____ nb_elem_std_ " << finl;
  os << nb_elem_std_ << finl;
  os << "____ volumes_entrelaces_ " << finl;
  volumes_entrelaces_.ecrit(os);
  os << "____ face_normales_ " << finl;
  face_normales_.ecrit(os);
  os << "____ facette_normales_ " << finl;
  facette_normales_.ecrit(os);
  os << "____ nb_faces_std_ " << finl;
  os << nb_faces_std_ << finl;
  os << "____ rang_elem_non_std_ " << finl;
  rang_elem_non_std_.ecrit(os);
  return os;
}
 
Sortie& Domaine_VEF::printOn(Sortie& os) const
{
  Domaine_VF::printOn(os);
 
  os << h_carre << finl;
  os << type_elem_.valeur() << finl;
  os << nb_elem_std_ << finl;
  os << volumes_entrelaces_ << finl;
  os << face_normales_ << finl;
  os << facette_normales_ << finl;
  os << xp_ << finl;
  os << xv_ << finl;
  os << nb_faces_std_ << finl;
  os << rang_elem_non_std_ << finl;
  return os;
}
 
Entree& Domaine_VEF::readOn(Entree& is)
{
  Domaine_VF::readOn(is);
  is >> h_carre;
 
  /* read type_elem */
  {
    Nom type;
    is >> type;
    if (type == "Tri_VEF")
      type_elem_ = Tri_VEF();
    else if (type == "Tetra_VEF")
      type_elem_ = Tetra_VEF();
    else if (type == "Quadri_VEF")
      type_elem_ = Quadri_VEF();
    else if (type == "Hexa_VEF")
      type_elem_ = Hexa_VEF();
    else
      {
        Cerr << type << " is not an Elem_VEF!" << finl;
        Process::exit();
      }
  }
 
  is >> nb_elem_std_;
  is >> volumes_entrelaces_;
  is >> face_normales_;
  is >> facette_normales_;
  is >> xp_;
  is >> xv_;
  is >> nb_faces_std_;
  is >> rang_elem_non_std_;
  return is;
}
 
void exemple_champ_non_homogene(const Domaine_VEF& domaine_VEF, DoubleTab& tab)
{
  const DoubleTab& xp = domaine_VEF.xp();
  const Domaine& domaine=domaine_VEF.domaine();
  const DoubleTab& coord=domaine.coord_sommets();
  const DoubleTab& xa=domaine_VEF.xa();
  const ArrOfInt& renum_arete_perio=domaine_VEF.get_renum_arete_perio();
  // Check array xa of edge coordinates
  if (xa.size_array()) Debog::verifier("xa=",xa);
  int nb_elem=domaine.nb_elem();
  int nb_elem_tot=domaine.nb_elem_tot();
  int nb_som=domaine.nb_som();
  int nb_som_tot=domaine.nb_som_tot();
  int nb_aretes=domaine.nb_aretes();
  for (int I=0 ; I<nb_elem; I++)
    {
      tab(I)=(1.1+xp(I,0))*(1.1+2*xp(I,1));
      if (Objet_U::dimension==3) tab(I)*=(1.1+3*xp(I,2));
    }
  for (int I=0; I<nb_som; I++)
    {
      tab(nb_elem_tot+I)=(1.1+coord(I,0))*(1.1+2*coord(I,1));
      if (Objet_U::dimension==3) tab(nb_elem_tot+I)*=(1.1+3*coord(I,2));
      // Apply periodicity:
      tab(nb_elem_tot+I)=tab(nb_elem_tot+domaine.get_renum_som_perio(I));
    }
 
#ifndef NDEBUG
  const IntVect& ok_arete = domaine_VEF.get_ok_arete();
#endif
  for (int I=0; I<nb_aretes; I++)
    {
      tab(nb_elem_tot+nb_som_tot+I)=(1.1+xa(I,0))*(1.1+2*xa(I,1))*(1.1+3*xa(I,2));
      // Apply periodicity:
      tab(nb_elem_tot+nb_som_tot+I)=tab[nb_elem_tot+nb_som_tot+renum_arete_perio[I]];
      // Check ok_arete in passing
      assert(ok_arete(I)==ok_arete(renum_arete_perio[I]));
    }
  tab.echange_espace_virtuel();
}
 
void Domaine_VEF::swap(int fac1, int fac2, int nb_som_faces)
{
 
}
 
void Domaine_VEF::discretiser()
{
  Domaine& domaine_geom = domaine();
  typer_elem(domaine_geom);
  Elem_geom_base& elem_geom = domaine_geom.type_elem().valeur();
  Domaine_VF::discretiser();
 
  // Correction of the face_voisins array:
  //  After Domaine_VF::discretiser(), neighbor elements 0 and 1 of a face
  //  are the same on all processors that own the face.
  //  If the face is virtual and one of its two neighbor elements is
  //  unknown (not in the ghost layer), the neighbor element
  //  is -1. This can be neighbor 0 or neighbor 1.
  //  Correct virtual faces so that if a neighbor element is unknown,
  //  it is neighbor1. Neighbor0 is therefore always valid.
  {
    IntTab& face_vois = face_voisins();
    const int debut = nb_faces();
    const int fin = nb_faces_tot();
    for (int i = debut; i < fin; i++)
      {
        if (face_voisins(i, 0) == -1)
          {
            face_vois(i, 0) = face_vois(i, 1);
            face_vois(i, 1) = -1;
          }
      }
  }
 
  // Check coherence between the geometric element and
  //the discretisation element
 
  if (sub_type(Tri_VEF, type_elem_.valeur()))
    {
      if (!sub_type(Triangle, elem_geom))
        {
          Cerr << " The geometric element type " << elem_geom.que_suis_je() << " is incorrect" << finl;
          Cerr << " Only the Triangle type is compatible with VEF discretization in 2D" << finl;
          Cerr << " The domain must be triangulated when using the internal mesher";
          Cerr << " using the instruction: Trianguler nom_dom" << finl;
          exit();
        }
    }
  else if (sub_type(Tetra_VEF, type_elem_.valeur()))
    {
      if (!sub_type(Tetraedre, elem_geom))
        {
          Cerr << " The geometric element type " << elem_geom.que_suis_je() << " is incorrect" << finl;
          Cerr << " Only the Tetraedre type is compatible with VEF discretization in 3D" << finl;
          Cerr << " The domain must be tetrahedralized when using the internal mesher";
          Cerr << " using the instruction: Tetraedriser nom_dom" << finl;
          exit();
        }
    }
 
  else if (sub_type(Quadri_VEF, type_elem_.valeur()))
    {
 
      if (!sub_type(Quadrangle_VEF, elem_geom))
        {
          Cerr << " The geometric element type " << elem_geom.que_suis_je() << " is incorrect" << finl;
          exit();
        }
    }
  else if (sub_type(Hexa_VEF, type_elem_.valeur()))
    {
 
      if (!sub_type(Hexaedre_VEF, elem_geom))
        {
          Cerr << " The geometric element type " << elem_geom.que_suis_je() << " is incorrect" << finl;
          exit();
        }
    }
 
  // 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.
  {
    const int n = nb_faces();
    face_normales_.resize(n, dimension);
    creer_tableau_faces(face_normales_, RESIZE_OPTIONS::NOCOPY_NOINIT);
    const IntTab& face_som = face_sommets();
    const IntTab& face_vois = face_voisins();
    const IntTab& elem_face = elem_faces();
    type_elem_->creer_face_normales(face_normales_, face_som, face_vois, elem_face, domaine_geom);
  }
 
  // Compute facette_normales_
  type_elem_->creer_facette_normales(*this, rang_elem_non_std());
 
  calculer_volumes_entrelaces();
  Cerr << "Informations of the Domaine VEF of the domain " << domaine().le_nom() << " : " << finl;
 
  calculer_h_carre();
 
  domaine().creer_tableau_sommets(volumes_som_, RESIZE_OPTIONS::NOCOPY_NOINIT);
 
  double coeff=1./3.;
  if (dimension==3)
    coeff=1./4.;
 
  const IntTab& elements = domaine().les_elems();
  const int nb_som_elem = elements.dimension(1);
  // Loop over all elements because contributions are added to joint vertices:
  const int n = nb_elem_tot();
  // The virtual space of volumes is needed
  const DoubleVect& volume_elem = volumes();
  assert_espace_virtuel_vect(volume_elem);
 
  // Zero out the entire array because += is done on virtual items
  // (otherwise uninitialized memory is accessed)
  operator_egal(volumes_som_, 0., VECT_ALL_ITEMS);
  for(int k=0; k<n; k++)
    {
      double volume = coeff * volume_elem(k);
      for(int isom=0; isom<nb_som_elem; isom++)
        {
          int som = elements(k, isom);
          volumes_som_(som)+=volume;
        }
    }
  volumes_som_.echange_espace_virtuel();
}
 
void Domaine_VEF::discretiser_suite(const VEF_discretisation& discr)
{
  // Retrieve discretisation parameters
  alphaE = discr.get_alphaE();
  alphaS = discr.get_alphaS();
  alphaA = discr.get_alphaA();
  alphaRT = discr.get_alphaRT();
  P1Bulle = discr.get_P1Bulle();
  modif_div_face_dirichlet= discr.get_modif_div_face_dirichlet();
  cl_pression_sommet_faible = discr.get_cl_pression_sommet_faible();
 
 
  if (alphaA)
    discretiser_arete();
 
  // Build the descriptor for p1bulle arrays
  {
    MD_Vector_composite md_p1b;
    if (alphaE)
      {
        const MD_Vector& md = domaine().md_vector_elements();
        md_p1b.add_part(md, 0, "P0");
      }
    if (alphaS)
      {
        const MD_Vector& md = domaine().md_vector_sommets();
        md_p1b.add_part(md, 0, "P1");
      }
    if (alphaA)
      {
        const MD_Vector& md = md_vector_aretes();
        md_p1b.add_part(md, 0, "Pa");
      }
    md_vector_p1b_.copy(md_p1b);
  }
  Cerr << "the Domaine_VEF has been filled successfully" << finl;
}
 
/*
  static int arete(const IntTab& som_aretes,
  const IntTab& aretes_som,
  int S1, int S2, const Domaine& dom)
  {
  //Cout << "S1,S2 = " << S1 << " " << S2 << finl;
  for(int i=0; i<nb_max_aretes_som; i++)
  {
  const int& a=som_aretes(S1,i);
  if( (dom.get_renum_som_perio(aretes_som(a,0))==S2)
  ||
  (dom.get_renum_som_perio(aretes_som(a,1))==S2) )
  return a;
  }
  {
  Cerr << "arete pas Trouvee!!" << finl;
  Cout << "S1,S2 = " << S1 << " " << S2 << finl;
  Cout << "som_aretes " << som_aretes << finl;
  Cout << "aretes_som " << aretes_som << finl;
  Process::exit();
  }
  return -1;
  } */
 
static int next(int S,
                const ArrOfInt& contenu)
{
  int nb_som=contenu.size_array();
  int i=S+1;
  while(i<nb_som)
    if(contenu[i]==1)
      return i;
    else i++;
  i=0;
  while(i<S)
    if(contenu[i]==1)
      return i;
    else i++;
  return -1;
}
 
void Domaine_VEF::discretiser_arete()
{
  const Domaine& dom = domaine();
 
  // Create real edges (geometric information built and stored in the domain)
  domaine().creer_aretes();
  md_vector_aretes_ = domaine().aretes_som().get_md_vector();
 
  // Compute edge gravity centers xa_ stored in Domaine_VF
  const IntTab& aretes_som = domaine().aretes_som();
  const int nb_aretes = aretes_som.dimension(0);
  const DoubleTab& coord = dom.les_sommets();
  const int dim = coord.dimension(1);
  xa_.resize(0, dim);
  creer_tableau_aretes(xa_, RESIZE_OPTIONS::NOCOPY_NOINIT);
  for (int i = 0; i < nb_aretes; i++)
    {
      const int s0 = aretes_som(i, 0);
      const int s1 = aretes_som(i, 1);
      for (int j = 0; j < dim; j++)
        xa_(i, j) = (coord(s0, j) + coord(s1, j)) * 0.5;
    }
  xa_.echange_espace_virtuel();
 
  const IntTab& elem_aretes = domaine().elem_aretes();
 
  // Compute edge volumes
  // Create an array initialised to zero:
  creer_tableau_aretes(volumes_aretes);
 
  const int nbr_elem = domaine().nb_elem();
  const int nb_aretes_elem = elem_aretes.dimension(1);
  // factor 6 for edge volume computation, is this correct for elements other than tetra?
  assert(nb_aretes_elem == 6);
  // Alternative approach: compute contributions to edges from real elements,
  // then sum contributions of shared edges:
  for (int elem = 0; elem < nbr_elem; elem++)
    {
      double vol = volumes(elem) / 6.0;
      for(int j = 0; j < nb_aretes_elem; j++)
        {
          int arete = elem_aretes(elem, j);
          volumes_aretes[arete] += vol;
        }
    }
  // Sum joint edge contributions and exchange virtual spaces
  MD_Vector_tools::echange_espace_virtuel(volumes_aretes, MD_Vector_tools::EV_SOMME_ECHANGE);
}
 
void Domaine_VEF::construire_ok_arete()
{
  Cerr << "Build array ok_arete..." << finl;
  const Domaine& dom=domaine();
  const IntTab& aretes_som=domaine().aretes_som();
  const int nb_som_reel=nb_som();
 
  // Vertex-to-edge connectivity (for each vertex, list of adjacent edges)
  // B.M.: replacing IntTab(n, 64) with a Static_Int_List, reusing
  //  the connectivity computation method... more memory efficient!
  // som_aretes contains a modified connectivity for periodic vertices
  // (edges are always attached to get_renum_som_perio() vertex,
  //  opposite periodic edges are replaced by a fictitious edge [nb_som_reel,nb_som_reel],
  //  and opposite periodic vertices are not attached to any edge)
  const int nb_aretes = aretes_som.dimension(0);
  Static_Int_Lists som_aretes;
  {
    // Create an edge array where vertices are replaced by their periodic counterpart
    //  and without the opposite periodic edges
    IntTab aretes_som2;
    aretes_som2.copy(aretes_som, RESIZE_OPTIONS::NOCOPY_NOINIT);
 
    for (int i = 0; i < nb_aretes; i++)
      {
        if (renum_arete_perio[i] == i)
          {
            aretes_som2(i, 0) = dom.get_renum_som_perio(aretes_som(i,0));
            aretes_som2(i, 1) = dom.get_renum_som_perio(aretes_som(i,1));
          }
        else
          {
            // opposite periodic edge, not retained; create a fictitious edge for
            // construire_connectivite_som_elem: vertex with index nb_som_reel
            // will be connected to all suppressed periodic edges
            aretes_som2(i, 0) = nb_som_reel;
            aretes_som2(i, 1) = nb_som_reel;
          }
      }
    // Add one extra vertex (the fictitious vertex connected to suppressed edges)
    construire_connectivite_som_elem(nb_som_reel+1, aretes_som2, som_aretes, 0 /* do not include virtual items */);
  }
 
  // Create and initialise array contenu, initialised to zero by default
  ArrOfInt contenu(nb_som_reel);
 
  // Parallel initialisation: contenu is set to 2
  // for shared vertices received from another processor
  {
    ArrOfBit flags;
    dom.les_sommets().get_md_vector()->get_sequential_items_flags(flags);
    for (int i = 0; i < nb_som_reel; i++)
      {
        if (!flags[i])
          {
            // This vertex is received from another processor
            const int i2 = dom.get_renum_som_perio(i);
            contenu[i] = 2;
            if (i2 != i)
              contenu[i2] = 2;
          }
      }
  }
 
  // Estimate the number of superfluous edges to find on the domain
  // by counting real non-periodic vertices whose contenu equals 0
  int nombre_aretes_superflues_prevues_sur_le_dom=0;
  for (int i=0; i<nb_som_reel; i++)
    if (i==dom.get_renum_som_perio(i) && contenu[i]==0)
      nombre_aretes_superflues_prevues_sur_le_dom++;
 
  ok_arete = -1;
 
  // Loop while there are vertices with contenu == 0
  while (min_array(contenu)==0)
    {
      // Find the first vertex with contenu == 0 by looping over edges
      int Aroot=-1,Sroot=-1;
      int S0,S1;
      do
        {
          while (ok_arete(++Aroot)==0) {};
          S0=dom.get_renum_som_perio(aretes_som(Aroot,0));
          S1=dom.get_renum_som_perio(aretes_som(Aroot,1));
        }
      while (Aroot<nb_aretes && contenu[S0]!=0 && contenu[S1]!=0);
 
      assert(Aroot<nb_aretes);
 
      Aroot=renum_arete_perio[Aroot];
 
      ok_arete(Aroot)=0; // lockee
      if (!contenu[S0])
        {
          contenu[S0]=1;
          Sroot=S0;
        }
      else if (!contenu[S1])
        {
          contenu[S1]=1;
          Sroot=S1;
        }
 
      // Loop over vertices
      do
        {
          const int nb_aretes_voisines = som_aretes.get_list_size(Sroot);
          for(int i=0; i<nb_aretes_voisines; i++)
            {
              int A=renum_arete_perio[som_aretes(Sroot,i)];
              if(ok_arete(A)==-1)
                {
                  int S=dom.get_renum_som_perio(aretes_som(A,0));
                  if(S==Sroot)
                    {
                      S=dom.get_renum_som_perio(aretes_som(A,1));
                    }
                  if(!contenu[S])
                    {
                      ok_arete(A)=0; // lockee
                      contenu[S]=1;
                    }
                  else
                    {
                      ok_arete(A)=1; // libre
                    }
                }
            }
          contenu[Sroot]=2;
        }
      while((Sroot=next(Sroot, contenu))!=-1);
 
      // Correction of ok_arete and contenu arrays for periodicity
      for(int i=0; i<nb_aretes; i++)
        ok_arete(i)=ok_arete[renum_arete_perio[i]];
      for(int i=0; i<nb_som_reel; i++)
        contenu[i]=contenu[dom.get_renum_som_perio(i)];
    }
 
  // Update virtual parts of the ok_arete array
  ok_arete.echange_espace_virtuel();
 
  // Check superfluous edges
  verifie_ok_arete(nombre_aretes_superflues_prevues_sur_le_dom);
 
  // Write superfluous edges to file nom_du_cas.ok_arete to be read next time
  Nom fichier(nom_du_cas());
  fichier+="_";
  fichier+=domaine().le_nom()+".ok_arete";
 
  Cerr << "Writing file " << fichier << finl;
  EcrFicPartageBin fic_ok_arete_;
  if (!fic_ok_arete_.ouvrir(fichier,ios::out))
    {
      Cerr << "Error: cannot open file " << fichier << " for writing." << finl;
      exit();
    }
 
  ArrOfBit marqueurs_aretes;
  const MD_Vector& md_aretes = md_vector_aretes();
  md_aretes->get_sequential_items_flags(marqueurs_aretes);
  const trustIdType nb_aretes_seq = md_aretes->nb_items_seq_tot();
 
  if (Process::je_suis_maitre())
    fic_ok_arete_ << nb_aretes_seq << finl;
 
  // Each processor writes its non-shared edges:
  const int n = marqueurs_aretes.size_array();
  for (int i = 0; i < n; i++)
    {
      if (marqueurs_aretes[i])
        fic_ok_arete_ << xa_(i,0) << tspace << xa_(i,1) << tspace << xa_(i,2) << tspace << ok_arete(i) << finl;
    }
  fic_ok_arete_.syncfile();
  fic_ok_arete_.close();
}
 
void Domaine_VEF::construire_renum_arete_perio(const Conds_lim& conds_lim)
{
  Cerr << "Build array renum_arete_perio..." << finl;
  const IntTab& aretes_som=domaine().aretes_som();
  const int nb_aretes_tot = static_cast<int>(domaine().nb_aretes_tot()); // domain is already discretised, so already split, so we're just working with a small part
  const Domaine& dom=domaine();
 
  // Initialise renum_arete_perio
  renum_arete_perio.resize_array(nb_aretes_tot);
  for (int i=0; i<nb_aretes_tot; i++)
    renum_arete_perio[i]=i;
 
  const IntTab& elem_aretes=domaine().elem_aretes();
  ArrOfInt aretes1(6);
  ArrOfInt aretes2(6);
  // First step: make all periodic edges linked together point to the same edge
  for (auto& itr : conds_lim)
    {
      //for cl
 
      const Cond_lim_base& cl = itr.valeur();
      if (sub_type(Periodique, cl))
        {
          //if Perio
          const Periodique& la_cl_perio = ref_cast(Periodique,cl);
          const Front_VF& le_bord = ref_cast(Front_VF,cl.frontiere_dis());
 
          int nf_bord_tot = le_bord.nb_faces_tot();
          IntVect fait(nf_bord_tot);
          fait = 0;
          for(int ind_face=0; ind_face<nf_bord_tot; ind_face++)
            if(!fait(ind_face))
              {
                int face=le_bord.num_face(ind_face);
                int ind_faassociee=la_cl_perio.face_associee(ind_face);
                int face_assciee=le_bord.num_face(ind_faassociee);
                fait(ind_faassociee)=fait(ind_face)=1;
 
                int elem1 = face_voisins_(face,0);
                int elem2 = face_voisins_(face,1);
 
                for(int j=0; j<6; j++)
                  {
                    aretes1[j]=elem_aretes(elem1,j);
                    aretes2[j]=elem_aretes(elem2,j);
                  }
 
                for(int j1=0; j1<6; j1++)
                  {
                    int ar1=aretes1[j1];
                    int& ar1_perio = renum_arete_perio[ar1];
                    // Check that the edge belongs to the face (ok==2)
                    int som11=aretes_som(ar1, 0);
                    int som12=aretes_som(ar1, 1);
                    int ok=0;
                    int nbf = domaine().type_elem()->nb_som_face();
                    const IntTab& sommet = face_sommets();
                    for (int k=0; k<nbf; k++)
                      if (sommet(face,k)==som11 || sommet(face,k)==som12) ok++;
                    assert(ok>0);
                    if (ok==2)
                      {
                        int s11=dom.get_renum_som_perio(som11);
                        int s12=dom.get_renum_som_perio(som12);
                        for(int j2=0; j2<6; j2++)
                          {
                            int ar2=aretes2[j2];
                            int& ar2_perio = renum_arete_perio[ar2];
                            int som21=aretes_som(ar2, 0);
                            int som22=aretes_som(ar2, 1);
                            ok=0;
                            // Check that the edge belongs to face_assciee (ok==2)
                            for (int k=0; k<nbf; k++)
                              if (sommet(face_assciee,k)==som21 || sommet(face_assciee,k)==som22) ok++;
                            assert(ok>0);
                            if (ok==2)
                              {
                                int s21=dom.get_renum_som_perio(som21);
                                int s22=dom.get_renum_som_perio(som22);
                                assert(ar1!=ar2);
                                //The edges are therefore periodic if we enter the "if"
                                if ( ( (s21==s11)||(s22==s11) ) && ( (s22==s12)||(s21==s12) ) )
                                  {
                                    // Criterion I: choose the periodic edge based on a geometric criterion
                                    int dir_perio = la_cl_perio.direction_periodicite();
                                    int arete_perio = (xa(ar1_perio,dir_perio)<=xa(ar2_perio,dir_perio)) ? ar1_perio : ar2_perio;
                                    /*
                                    // Criterion II: choose the periodic edge based on periodic vertices
                                    int arete_perio;
                                    if (som11!=s11 && som12!=s12) // Vertices of edge 1 are periodic
                                    arete_perio = ar2_perio;
                                    else if (som21!=s21 && som22!=s22) // Vertices of edge 2 are periodic
                                    arete_perio = ar1_perio;
                                    else
                                    {
                                    Cerr << "Cas non prevu." << finl;
                                    exit();
                                    } */
                                    ar2_perio = arete_perio;
                                    ar1_perio = arete_perio;
                                    volumes_aretes(ar2)+=volumes_aretes(ar1);
                                    volumes_aretes(ar1)+=volumes_aretes(ar2);
                                  }//end if "s21==s11"
                              }
                          }//end for "j2"
                      }
                  }//end if "renum_arete_perio(aretes1)" and end for "j1"
              }//end if "fait"
        }//end if Perio
    }// end for cl
 
  // Second step: make all periodic edges linked together point to the same edge
  for (int i = 0; i < nb_aretes_tot; i++)
    {
      int j = renum_arete_perio[i];
      // Traverse linked edges for this chain:
      while (j != renum_arete_perio[j])
        j = renum_arete_perio[j];
      renum_arete_perio[i] = j;
    }
 
  if (Debog::active())
    {
      IntVect tmp;
      creer_tableau_aretes(tmp, RESIZE_OPTIONS::NOCOPY_NOINIT);
      const int n = tmp.size_array();
      for (int i=0; i<n; i++)
        tmp[i] = renum_arete_perio[i];
      Debog::verifier_indices_items("renum_arete_perio",tmp.get_md_vector(),tmp);
    }
}
 
void Domaine_VEF::verifie_ok_arete(int nombre_aretes_superflues_prevues_sur_le_dom) const
{
  Cerr << "Check array ok_arete..." << finl;
  // Algorithm to verify the superfluous edge array
  // ok_arete(i)==0 : edge i is superfluous
  // ok_arete(i)==2 : edge i is necessary
  // ...
  // contenu must contain only 2, meaning all vertices have been analysed.
  const Domaine& dom=domaine();
  const int nb_som_reel=nb_som();
  const IntTab& aretes_som=domaine().aretes_som();
  int nb_aretes_pour_verbose=60; // for verbose output
  ArrOfInt sommet_relie_arete_superflue(nb_som_reel);
  sommet_relie_arete_superflue=0;
  for (int i=0; i<nb_som_reel; i++) // If i is a periodic vertex:
    if (dom.get_renum_som_perio(i)!=i) // it suffices to verify on vertex dom.get_renum_som_perio(i)
      sommet_relie_arete_superflue[i]=1; // so we skip i
 
  double nombre_aretes_reelles_superflues=0;
  // Traverse all edges but only consider real vertices
  int nb_aretes_tot=domaine().nb_aretes_tot();
  int nb_aretes_reelles=domaine().nb_aretes();
  for (int i=0; i<nb_aretes_tot; i++)
    {
      if (!ok_arete(i)) // arete superflue
        {
          int S0=aretes_som(i,0);
          if (S0<nb_som_reel)
            {
              sommet_relie_arete_superflue[S0]=1;
              sommet_relie_arete_superflue[dom.get_renum_som_perio(S0)]=1;
            }
          int S1=aretes_som(i,1);
          if (S1<nb_som_reel)
            {
              sommet_relie_arete_superflue[S1]=1;
              sommet_relie_arete_superflue[dom.get_renum_som_perio(S1)]=1;
            }
          if (renum_arete_perio[i]==i) // Count periodic edges only once
            {
              if (nb_aretes_tot<nb_aretes_pour_verbose)
                {
                  if (renum_arete_perio[i]==i) Cerr << "[" << Process::me() << "] Arete " << i << " superflue non perio: " << S0 << " " << S1 << finl;
                  else Cerr << "[" << Process::me() << "] Arete " << i << " superflue perio: " << S0 << " " << S1 << " Periodique avec " << renum_arete_perio[i] << finl;
                }
              int aretes_superflues_communes=1;
              // If the superfluous edge is shared, account for it
              for (int j=0; j<dom.faces_joint().size(); j++)
                {
                  int nb_aretes_sur_le_joint = dom.faces_joint()(j).joint_item(JOINT_ITEM::ARETE).items_communs().size_array();
                  for (int k=0; k<nb_aretes_sur_le_joint; k++)
                    if (dom.faces_joint()(j).joint_item(JOINT_ITEM::ARETE).items_communs()[k]==i) aretes_superflues_communes++;
                }
              // Count real superfluous edges
              if (i<nb_aretes_reelles)
                nombre_aretes_reelles_superflues+=1./aretes_superflues_communes;
            }
        }
    }
 
  if (nb_som_reel > 0 && min_array(sommet_relie_arete_superflue)==0)
    {
      Cerr << finl << "[" << Process::me() << "] There is at least one vertex not linked to a redundant edge:" << finl;
      for (int i=0; i<sommet_relie_arete_superflue.size_array(); i++)
        if (sommet_relie_arete_superflue[i]==0) Cerr << "Vertex " << i << finl;
      Process::exit();
    }
  Cerr << "[" << Process::me() << "] Verification that each non-periodic vertex is linked to at least one redundant edge: OK!" << finl;
  // Count the total number of vertices accounting for periodicity and shared vertices
  double nb_sommets_non_periodiques=0;
  for (int i=0; i<nb_som_reel; i++)
    if (dom.get_renum_som_perio(i)==i) // Non-periodic vertex
      {
        int sommets_communs=1;
        // Account for shared vertices
        for (int j=0; j<dom.faces_joint().size(); j++)
          for (int k=0; k<dom.faces_joint()(j).joint_item(JOINT_ITEM::SOMMET).items_communs().size_array(); k++)
            if (dom.faces_joint()(j).joint_item(JOINT_ITEM::SOMMET).items_communs()[k]==i) sommets_communs++;
        nb_sommets_non_periodiques+=1./sommets_communs;
      }
  double total_nombre_aretes_superflues = mp_sum(nombre_aretes_reelles_superflues);
  double somme_nombre_aretes_superflues_prevues_par_domaine = mp_sum_as_double(nombre_aretes_superflues_prevues_sur_le_dom);
  double total_nb_sommets_non_periodiques = mp_sum(nb_sommets_non_periodiques);
 
 
  // Cerr << "Nombre de sommets non periodiques           = " << total_nb_sommets_non_periodiques << finl;
  // Cerr << "Nombre de sommets periodiques               =  " << nb_som_reel-total_nb_sommets_non_periodiques << finl;
  //int nb_aretes_periodiques=0;
  //int nb_aretes_perio_superflues=0;
  for (int i=0; i<nb_aretes_tot; i++)
    {
      if (renum_arete_perio[i]!=i)
        {
          //nb_aretes_periodiques++;
          assert(ok_arete(i)==ok_arete(renum_arete_perio[i]));
          //if (!ok_arete(i)) nb_aretes_perio_superflues++;
        }
    }
  // Cerr << "Number of non-periodic edges             = " << nb_aretes_tot-nb_aretes_periodiques << finl;
  // Cerr << "Number of periodic edges                 =  " << nb_aretes_periodiques << finl;
 
  // Cerr << "Number of redundant edges                =  " << total_nombre_aretes_superflues << finl;
  // Cerr << "Number of periodic redundant edges       =  " << nb_aretes_perio_superflues << finl;
  // Cerr << "Number of non-redundant periodic edges   =  " << nb_aretes_periodiques-nb_aretes_perio_superflues << finl;
 
  if (Process::is_sequential()) // Limited to sequential because there is an issue with computing total_nombre_aretes_superflues (common items for edges are not yet built!)
    if (!est_egal(somme_nombre_aretes_superflues_prevues_par_domaine,total_nombre_aretes_superflues) && je_suis_maitre())
      {
        Cerr << "The sum of redundant edges predicted per domain is not equal to the number of redundant edges found on the domain." << finl;
        Cerr << somme_nombre_aretes_superflues_prevues_par_domaine << " != " << total_nombre_aretes_superflues << finl;
        Process::exit();
      }
  Cerr << "Total number of redundant edges: " << somme_nombre_aretes_superflues_prevues_par_domaine << finl;
  Cerr << "Total number of non-periodic vertices = " << total_nb_sommets_non_periodiques << finl;
  Cerr << "Verifying equality between total number of redundant edges and total number of vertices: ";
  // Verify that the number of redundant edges equals the number of vertices
  if (!est_egal(somme_nombre_aretes_superflues_prevues_par_domaine,total_nb_sommets_non_periodiques) && je_suis_maitre())
    {
      Cerr << "Mismatch. Failure of the redundant edge search algorithm." << finl;
      Process::exit();
    }
  Cerr << "OK!" << finl << "Redundant edge verification correct." << finl;
}
 
// Return value:
//  1: ok
//  0: file does not exist
//  -1: file exists but wrong number of edges
int Domaine_VEF::lecture_ok_arete()
{
  Nom fichier(nom_du_cas());
  fichier+="_";
  fichier+=domaine().le_nom()+".ok_arete";
 
  Cerr << "Trying to read file " << fichier << " (edges to remove from the set of degrees of freedom)" << finl;
  EFichierBin fic_ok_arete_;
  // Read ok_arete from file for sequential-parallel comparison
  if (!fic_ok_arete_.ouvrir(fichier,ios::out))
    {
      Cerr << "File " << fichier << " does not exist." << finl;
      return 0;
    }
 
  int n;
  fic_ok_arete_ >> n;
  if (n != md_vector_aretes()->nb_items_seq_tot())
    {
      Cerr << "File " << fichier << " is not compatible with the current mesh." << finl;
      return 0;
    }
 
  Octree_Double octree;
  octree.build_nodes(xa_, 1 /* include virtual nodes */);
  const double eps = precision_geom;
  ArrOfInt liste_aretes;
 
 
  IntVect marqueurs;
  creer_tableau_aretes(marqueurs); // initialised to zero by default
 
  for (int i = 0; i < n; i++)
    {
      double x, y, z;
      int flag;
      fic_ok_arete_ >> x >> y >> z >> flag;
      octree.search_elements_box(x-eps, y-eps, z-eps, x+eps, y+eps, z+eps, liste_aretes);
      octree.search_nodes_close_to(x, y, z, xa_, liste_aretes, eps);
      const int m = liste_aretes.size_array();
      if (m > 1)
        {
          Cerr << "File " << fichier << " is not compatible with the current mesh." << finl;
          return 0;
        }
      if (m == 1)
        {
          const int arete = liste_aretes[0];
          ok_arete[arete] = flag;
          marqueurs[arete]++;
        }
      else
        {
          // Normal in parallel: not all edges are available on all processors...
        }
    }
 
  // Each edge must have been found exactly once:
  if (max_array(marqueurs) > 1 || min_array(marqueurs) < 1)
    {
      Cerr << "File " << fichier << " is not compatible with the current mesh." << finl;
      return 0;
    }
  else
    Cerr << "File " << fichier << " is OK." << finl;
  fic_ok_arete_.close();
  assert_espace_virtuel_vect(ok_arete);
 
  return 1;
}
 
void Domaine_VEF::creer_tableau_p1bulle(Array_base& x, RESIZE_OPTIONS opt) const
{
  const MD_Vector& md = md_vector_p1b();
  MD_Vector_tools::creer_tableau_distribue(md, x, opt);
}
 
void Domaine_VEF::calculer_h_carre()
{
  const int nbe = nb_elem();
  // Compute h_carre
  h_carre = 1.e30;
  h_carre_.resize(nbe);
  // Compute face surfaces
  const int nb_faces_elem = domaine().nb_faces_elem();
  CDoubleArrView face_surfaces_v = face_surfaces().view_ro();
  CDoubleArrView volumes_v = volumes().view_ro();
  CIntTabView elem_faces_v = elem_faces().view_ro();
  DoubleArrView h_carre_v = h_carre_.view_rw();
  Kokkos::parallel_reduce(start_gpu_timer(__KERNEL_NAME__), range_1D(0, nbe),
                          KOKKOS_LAMBDA(const int num_elem, double& h_carre_local)
  {
    double surf_max = 0;
    for (int i = 0; i < nb_faces_elem; i++)
      {
        double surf = face_surfaces_v(elem_faces_v(num_elem, i));
        surf_max = (surf > surf_max) ? surf : surf_max;
      }
    double vol = volumes_v(num_elem) / surf_max;
    vol *= vol;
    h_carre_v(num_elem) = vol;
    if (vol < h_carre_local) h_carre_local = vol;
  }, Kokkos::Min<double>(h_carre));
  end_gpu_timer(__KERNEL_NAME__);
  h_carre = mp_min(h_carre);
  Cerr << "Lowest cell size h=(Volume/max(Surface))= " << sqrt(h_carre) << finl;
  double moyenne = 0.;
  Kokkos::parallel_reduce(start_gpu_timer(__KERNEL_NAME__), range_1D(0, nbe),
                          KOKKOS_LAMBDA(const int i, double& sum)
  {
    sum += h_carre_v(i);
  }, moyenne);
  end_gpu_timer(__KERNEL_NAME__);
  moyenne = mp_sum(moyenne);
  double h_carre_moyen = moyenne / mp_sum((double) nbe);
  Cerr << "Average cell size <h>= " << sqrt(h_carre_moyen) << finl;
  if (h_carre * 1e6 < h_carre_moyen)
    Cerr << "Warning, a cell size is more than 1000 times smaller than the average cell size of the domain. Check your mesh." << finl;
  Cerr << "==============================================" << finl;
}
 
void Domaine_VEF::calculer_volumes_entrelaces()
{
  //  Cerr << "face normals " << face_normales() << finl;
  // Compute interlaced volumes;
 
  // For a dynamic domain, the array may or may not already have the correct structure.
  // If not deformable, this function should not be called twice.
  assert(domaine().deformable() || !(volumes_entrelaces_.get_md_vector()));
  if (!(volumes_entrelaces_.get_md_vector() == md_vector_faces()))
    {
      volumes_entrelaces_.reset();
      creer_tableau_faces(volumes_entrelaces_, RESIZE_OPTIONS::NOCOPY_NOINIT);
    }
  int nb_faces_elem = domaine().nb_faces_elem();
  int num_face;
  int nmax = premiere_face_int();
  const double facteur = 1. / nb_faces_elem;
  for (num_face = 0; num_face < nmax; num_face++)
    {
      int elem1 = face_voisins_(num_face, 0);
      volumes_entrelaces_[num_face] = volumes(elem1) * facteur;
    }
  nmax = nb_faces();
  for (; num_face < nmax; num_face++)
    {
      int elem1 = face_voisins_(num_face, 0);
      int elem2 = face_voisins_(num_face, 1);
      volumes_entrelaces_[num_face] = (volumes(elem1) + volumes(elem2)) * facteur;
    }
  volumes_entrelaces_.echange_espace_virtuel();
}
 
void Domaine_VEF::modifier_pour_Cl(const Conds_lim& conds_lim)
{
  Journal() << "Domaine_VEF::Modifier_pour_Cl" << finl;
  for (auto &itr : conds_lim)
    {
      const Cond_lim_base& cl = itr.valeur();
      if (sub_type(Periodique, cl))
        {
          const Periodique& la_cl_period = ref_cast(Periodique, cl);
          int nb_faces_elem = domaine().nb_faces_elem();
          const Front_VF& la_front_dis = ref_cast(Front_VF, cl.frontiere_dis());
          // Modification of arrays face_voisins_, face_normales_, volumes_entrelaces_
          // Orient certain normals so that normals at periodic faces point
          // from face_voisins(face,0) to face_voisins(face,1),
          // consistent with internal faces
          ToDo_Kokkos("critical");
          for (int ind_face = 0; ind_face <  la_front_dis.nb_faces_tot(); ind_face++)
            {
              int face = la_front_dis.num_face(ind_face);
              if ((face_voisins_(face, 0) == -1) || (face_voisins_(face, 1) == -1))
                {
                  int faassociee = la_front_dis.num_face(la_cl_period.face_associee(ind_face));
                  int elem1 = face_voisins_(face, 0);
                  int elem2 = face_voisins_(faassociee, 0);
                  double vol = (volumes(elem1) + volumes(elem2)) / nb_faces_elem;
                  volumes_entrelaces_[face] = vol;
                  volumes_entrelaces_[faassociee] = vol;
                  face_voisins_(face, 1) = elem2;
                  face_voisins_(faassociee, 0) = elem1;
                  face_voisins_(faassociee, 1) = elem2;
                  double psc = 0;
                  for (int k = 0; k < dimension; k++)
                    psc += face_normales_(face, k) * (xv_(face, k) - xp_(face_voisins_(face, 0), k));
 
                  if (psc < 0)
                    for (int k = 0; k < dimension; k++)
                      face_normales_(face, k) *= -1;
 
                  for (int k = 0; k < dimension; k++)
                    face_normales_(faassociee, k) = face_normales_(face, k);
                }
            }
        }
    }
 
  // PQ: 10/10/05: periodic faces having double contribution,
  //               marquer_faces_double_contrib is called in this method
  //               to be able to use conds_lim.
  Domaine_VF::marquer_faces_double_contrib(conds_lim);
  // Build array num_fac_loc_
  construire_num_fac_loc();
 
  static DoubleVect* ptr=0;
  if(ptr!=&volumes_som_)
    {
      const Domaine& dom=domaine();
      const int ns = nb_som();
      CIntArrView renum_som_perio = dom.get_renum_som_perio().view_ro();
      DoubleArrView volumes_som = volumes_som_.view_rw();
      Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(0, ns), KOKKOS_LAMBDA(const int i)
      {
        int j = renum_som_perio(i);
        if (i != j)
          Kokkos::atomic_add(&volumes_som(j), volumes_som(i));
      });
      end_gpu_timer(__KERNEL_NAME__);
      Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(0, ns), KOKKOS_LAMBDA(const int i)
      {
        int j = renum_som_perio(i);
        if (i != j)
          volumes_som(i) = volumes_som(j);
      });
      end_gpu_timer(__KERNEL_NAME__);
      volumes_som_.echange_espace_virtuel();
      ptr = &volumes_som_;
    }
 
  // Check array renum_som_perio
  if (Debog::active())
    {
      IntVect tmp;
      const Domaine& dom = domaine();
      dom.creer_tableau_sommets(tmp, RESIZE_OPTIONS::NOCOPY_NOINIT);
      const int n = tmp.size_array();
      for (int i=0; i<n; i++)
        tmp[i] = dom.get_renum_som_perio(i);
      Debog::verifier_indices_items("renum_som_perio",tmp.get_md_vector(),tmp);
    }
 
  if(get_alphaA())
    {
      // Build renum_arete_perio
      construire_renum_arete_perio(conds_lim);
 
      // Create a distributed array for ok_arete
      creer_tableau_aretes(ok_arete, RESIZE_OPTIONS::NOCOPY_NOINIT);
 
      // If P1, find superfluous edges:
      if (get_alphaS())
        {
          // Try to read a .ok_arete file from a previous calculation
          // If it does not exist or is inconsistent with the mesh, rebuild ok_arete
          if (!lecture_ok_arete())
            construire_ok_arete();
        }
      else
        {
          // All edges are necessary if no P1 support
          ok_arete=1;
        }
      Debog::verifier_getref("ok_arete", ok_arete, ok_arete);
 
      //Debog::verifier("ok_arete (identical only if ok_arete re-read from the .ok_arete file):",ok_arete);
    }// end if 3D
}
 
void Domaine_VEF::typer_elem(Domaine& domaine_geom)
{
  const Elem_geom_base& elem_geom = domaine_geom.type_elem().valeur();
  if (sub_type(Rectangle, elem_geom))
    domaine_geom.typer("Quadrangle");
  else if (sub_type(Hexaedre, elem_geom))
    domaine_geom.typer("Hexaedre_VEF");
 
  const Nom& type_elem_geom = domaine_geom.type_elem()->que_suis_je();
 
  if (Motcle(type_elem_geom) != "Segment")
    {
      Nom type;
      if (type_elem_geom == "Triangle")
        type = "Tri_VEF";
      else if (type_elem_geom == "Tetraedre")
        type = "Tetra_VEF";
      else if (type_elem_geom == "Quadrangle")
        type = "Quadri_VEF";
      else if (type_elem_geom == "Hexaedre_VEF")
        type = "Hexa_VEF";
      else
        {
          Cerr << "typing problem in Elem_VEF::typer" << finl;
          Process::exit();
        }
      type_elem_.typer(type);
    }
}
 
DoubleTab& Domaine_VEF::vecteur_face_facette()
{
  // Build if size is zero
  // or if the mesh is deformable
  if (vecteur_face_facette_.size() == 0 || domaine().deformable())
    {
      // Size 8*n*4*3*2=192n
      const int nfa7 = type_elem().nb_facette();
      const int nb_poly_tot = nb_elem_tot();
      vecteur_face_facette_.resize(nb_poly_tot, nfa7, dimension, 2);
      const IntTab& KEL = type_elem().KEL();
      const IntTab& les_Polys = domaine().les_elems();
      const DoubleTab& coord = domaine().coord_sommets();
      const DoubleTab& xg = xp();
      int nb_som_facette=dimension;
      for (int poly = 0; poly < nb_poly_tot; poly++)
        {
          // Loop over facets of the non-standard polyhedron:
          for (int fa7 = 0; fa7 < nfa7; fa7++)
            {
              int num1 = elem_faces(poly, KEL(0, fa7));
              int num2 = elem_faces(poly, KEL(1, fa7));
 
              // Compute rx0 and rx1:
              for (int i = 0; i < dimension; i++)
                {
                  // Compute the i-th coordinate of the fa7 center
                  double coord_centre_fa7 = xg(poly, i);
                  for (int num_som_fa7 = 0; num_som_fa7 < nb_som_facette - 1; num_som_fa7++)
                    {
                      int isom_loc = KEL(num_som_fa7 + 2, fa7);
                      int isom_glob = les_Polys(poly, isom_loc);
                      coord_centre_fa7 += coord(isom_glob, i);
                    }
                  coord_centre_fa7 /= nb_som_facette;
                  // End compute i-th coordinate of fa7 center
                  vecteur_face_facette_(poly,fa7,i,0) = coord_centre_fa7 - xv_(num1,i);
                  vecteur_face_facette_(poly,fa7,i,1) = coord_centre_fa7 - xv_(num2,i);
                }
              // End of rx0 and rx1 computation
            }
        }
      Cerr << "Build of vecteur_face_facette() size:" << vecteur_face_facette_.size_array() << finl;
    }
  return vecteur_face_facette_;
}