Domaine_EF.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 <interface_INITGAUSS.h>
#include <interface_CALCULBIJ.h>
#include <Domaine_Cl_dis_base.h>
#include <Quadrangle_VEF.h>
#include <Domaine_Cl_EF.h>
#include <Equation_base.h>
#include <Hexaedre_VEF.h>
#include <Milieu_base.h>
#include <Domaine_EF.h>
#include <Periodique.h>
#include <Segment_EF.h>
#include <Rectangle.h>
#include <Tetraedre.h>
#include <Quadri_EF.h>
#include <Hexaedre.h>
#include <Triangle.h>
#include <Tetra_EF.h>
#include <EFichier.h>
#include <Segment.h>
#include <Hexa_EF.h>
#include <Domaine.h>
#include <Scatter.h>
#include <Tri_EF.h>
 
Implemente_instanciable(Domaine_EF, "Domaine_EF", Domaine_VF);
 
Sortie& Domaine_EF::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 << "____ 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_EF::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 << xp_ << finl;
  os << xv_ << finl;
  os << nb_faces_std_ << finl;
  os << rang_elem_non_std_ << finl;
  return os;
}
 
Entree& Domaine_EF::readOn(Entree& is)
{
  Domaine_VF::readOn(is);
  is >> h_carre;
 
  /* read type_elem */
  {
    Nom type;
    is >> type;
    if (type == "Tri_EF")
      type_elem_ = Tri_EF();
    else if (type == "Tetra_EF")
      type_elem_ = Tetra_EF();
    else if (type == "Quadri_EF")
      type_elem_ = Quadri_EF();
    else if (type == "Hexa_EF")
      type_elem_ = Hexa_EF();
    else
      {
        Cerr << type << " is not an Elem_EF !" << finl;
        Process::exit();
      }
  }
 
  is >> nb_elem_std_ ;
  is >> volumes_entrelaces_ ;
  is >> face_normales_ ;
  is >> xp_;
  is >> xv_;
  is >> nb_faces_std_ ;
  is >> rang_elem_non_std_ ;
  return is;
}
 
void Domaine_EF::calculer_IPhi(const Domaine_Cl_dis_base& zcl)
{
  int nbelem=domaine().nb_elem();
  int nb_som_elem=domaine().nb_som_elem();
  if (IPhi_.size() == 0) IPhi_.resize(nbelem,nb_som_elem);
  domaine().creer_tableau_elements(IPhi_);
  //  Scatter::creer_tableau_distribue(domaine().domaine(), JOINT_ITEM::ELEMENT, IPhi_);
  IPhi_thilde_=IPhi_;
  double rap=1./domaine().nb_som_elem();
 
  // computation of IPhi_ version valid for squares...
  for (int elem=0; elem<nbelem; elem++)
    {
      double vol=volumes(elem)*rap;
      for (int s=0; s<nb_som_elem; s++)
        {
          IPhi_(elem,s)=vol;
          IPhi_thilde_(elem,s)=vol*zcl.equation().milieu().porosite_elem(elem);
        }
    }
  IPhi_.echange_espace_virtuel();
  IPhi_thilde_.echange_espace_virtuel();
}
 
void Domaine_EF::calculer_volumes_sommets(const Domaine_Cl_dis_base& zcl)
{
  //  volumes_sommets_thilde_.resize(nb_som());
  //Scatter::creer_tableau_distribue(domaine().domaine(), JOINT_ITEM::SOMMET, volumes_sommets_thilde_);
//  domaine().creer_tableau_sommets(volumes_sommets_thilde_);
 
  calculer_IPhi(zcl);
 
  // computation of volumes_sommets_thilde
  const IntTab& les_elems=domaine().les_elems() ;
 
  DoubleVect volumes_som_prov(volumes_sommets_thilde_);
  for (int elem=0; elem<domaine().nb_elem_tot(); elem++)
    {
 
      for (int s=0; s<domaine().nb_som_elem(); s++)
        {
          double vol=IPhi_thilde_(elem,s);
          int s_glob=les_elems(elem,s);
          volumes_sommets_thilde_(s_glob)+=vol;
          volumes_som_prov(s_glob)+=IPhi_(elem,s);
        }
    }
 
//  volumes_sommets_thilde_.resize(nb_som());
  volumes_sommets_thilde_.echange_espace_virtuel();
  if (porosite_sommets_.size_array()==0)
    {
      // abort();
      porosite_sommets_=volumes_sommets_thilde_;
 
      for (int som=0; som<nb_som(); som++)
        {
          porosite_sommets_(som)=volumes_sommets_thilde_(som)/volumes_som_prov(som);
        }
      porosite_sommets_.echange_espace_virtuel();
    }
#if 0
  if (0)
    {
      porosite_sommets_.resize(0);
      DoubleVect volumes_som_prov(volumes_sommets_thilde_);
      DoubleVect volumes_som_prov_b(volumes_sommets_thilde_);
      volumes_som_prov=0;
      volumes_som_prov_b=0;
      for (int elem=0; elem<domaine().nb_elem_tot(); elem++)
        {
 
          for (int s=0; s<domaine().nb_som_elem(); s++)
            {
              double vol=IPhi_thilde_(elem,s);
              double v2=IPhi_(elem,s);
              int s_glob=les_elems(elem,s);
              volumes_som_prov(s_glob)+=v2/vol;
              volumes_som_prov_b(s_glob)+=1;
            }
        }
      Cerr<<mp_max_vect(porosite_elem_)<<" "<<mp_min_vect( porosite_elem_)<<finl;
      Cerr<<mp_max_vect(volumes_som_prov)<<" "<<mp_min_vect( volumes_som_prov)<<finl;
      Cerr<<mp_max_vect(volumes_som_prov_b)<<" "<<mp_min_vect( volumes_som_prov_b)<<finl;
 
      if (porosite_sommets_.size_array()==0)
        {
          // abort();
          porosite_sommets_=volumes_sommets_thilde_;
 
          for (int som=0; som<nb_som(); som++)
            {
              porosite_sommets_(som)=1./volumes_som_prov(som)*volumes_som_prov_b(som);
 
            }
          porosite_sommets_.echange_espace_virtuel();
        }
      Cerr<<mp_max_vect( porosite_sommets_)<<" "<<mp_min_vect( porosite_sommets_)<<finl;
      for (int elem=0; elem<domaine().nb_elem_tot(); elem++)
        for (int s=0; s<domaine().nb_som_elem(); s++)
          IPhi_thilde_(elem,s)= IPhi_thilde_(elem,s)*porosite_sommets_(les_elems(elem,s));
    }
#endif
 
  volumes_thilde_=volumes_;
  for (int elem=0; elem<domaine().nb_elem_tot(); elem++)
    {
      volumes_thilde_(elem)=volumes_(elem)*zcl.equation().milieu().porosite_elem(elem);
    }
}
 
void Domaine_EF::swap(int fac1, int fac2, int nb_som_faces )
{
 
}
 
void Domaine_EF::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();
  Nom type;
 
  if (type_elem_geom == "Triangle")
    type = "Tri_EF";
  else if (type_elem_geom == "Tetraedre")
    type = "Tetra_EF";
  else if (type_elem_geom == "Quadrangle")
    type = "Quadri_EF";
  else if (type_elem_geom == "Hexaedre_VEF")
    type = "Hexa_EF";
  else if (type_elem_geom == "Segment")
    type = "Segment_EF";
  else if (type_elem_geom == "Segment_axi")
    type = "Segment_EF_axi";
  else if (type_elem_geom == "Point")
    type = "Point_EF";
  else
    {
      Cerr << "typing problem in Domaine_EF::typer_elem => geometric type : " << type_elem_geom << finl;
      Process::exit();
    }
  type_elem_.typer(type);
}
 
void Domaine_EF::verifie_compatibilite_domaine()
{
  if (domaine().axi1d())
    {
      Cerr << "*****************************************************************************" << finl;
      Cerr << " Error in " << que_suis_je() << " : the type of domain " << domaine().que_suis_je();
      Cerr << " is not compatible" << finl;
      Cerr << " with the discretisation EF. " << finl;
      Cerr << " Please use the discretization EF_axi or define a domain of type Domaine." << finl;
      Cerr << "*****************************************************************************" << finl;
      Process::exit();
    }
}
 
void Domaine_EF::discretiser()
{
  verifie_compatibilite_domaine();
 
  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 the two neighbor elements is not
  //  known (it is not within the joint thickness), that neighbor element
  //  is -1. It can be either neighbor 0 or neighbor 1.
  //  Virtual faces are corrected so that, if a neighbor element is not
  //  known, it is placed as neighbor 1. Neighbor 0 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;
          }
      }
  }
 
  // Verification of coherence between the geometric element and
  // the discretisation element
 
 
  if (sub_type(Segment_EF,type_elem_.valeur()))
    {
      if (!sub_type(Segment,elem_geom))
        {
          Cerr << " The type of geometric element " << elem_geom.que_suis_je() << " is incorrect" << finl;
          exit();
        }
    }
  else if (sub_type(Tri_EF,type_elem_.valeur()))
    {
      if (!sub_type(Triangle,elem_geom))
        {
          Cerr << " The type of geometric element " <<
               elem_geom.que_suis_je() << " is incorrect" << finl;
          Cerr << " Only the Triangle type is compatible with the EF discretisation in dimension 2" << finl;
          Cerr << " You must triangulate the domain when using the internal mesher" ;
          Cerr << " using the instruction: Trianguler nom_dom" << finl;
          exit();
        }
    }
  else if (sub_type(Tetra_EF,type_elem_.valeur()))
    {
      if (!sub_type(Tetraedre,elem_geom))
        {
          Cerr << " The type of geometric element " <<
               elem_geom.que_suis_je() << " is incorrect" << finl;
          Cerr << " Only the Tetraedre type is compatible with the EF discretisation in dimension 3" << finl;
          Cerr << " You must tetrahedralise the domain when using the internal mesher";
          Cerr << " using the instruction: Tetraedriser nom_dom" << finl;
          exit();
        }
    }
 
  else if (sub_type(Quadri_EF,type_elem_.valeur()))
    {
 
      if (!sub_type(Quadrangle_VEF,elem_geom))
        {
          Cerr << " The type of geometric element " << elem_geom.que_suis_je() << " is incorrect" << finl;
          exit();
        }
    }
  else if (sub_type(Hexa_EF,type_elem_.valeur()))
    {
 
      if (!sub_type(Hexaedre_VEF,elem_geom))
        {
          Cerr << " The type of geometric element " << elem_geom.que_suis_je() << " is incorrect" << finl;
          exit();
        }
    }
 
  int num_face;
 
 
  // 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);
    // const Domaine & dom = domaine();
    //    Scatter::creer_tableau_distribue(dom, JOINT_ITEM::FACE, face_normales_);
    creer_tableau_faces(face_normales_);
    const IntTab& face_som = face_sommets();
    const IntTab& face_vois = face_voisins();
    const IntTab& elem_face = elem_faces();
    const int n_tot = nb_faces_tot();
    for (num_face = 0; num_face < n_tot; num_face++)
      {
        type_elem_->normale(num_face,
                            face_normales_,
                            face_som, face_vois, elem_face, domaine_geom);
      }
  }
 
  domaine().creer_tableau_sommets(volumes_sommets_thilde_);
}
 
void Domaine_EF::calculer_Bij_gen(DoubleTab& bij)
{
  const IntTab& les_elems=domaine().les_elems() ;
  int nbsom_face=nb_som_face();
  int nbelem=nb_elem_tot();
  int nbsom_elem=domaine().nb_som_elem();
  const IntTab& elemfaces=elem_faces_;
  int nbface_elem=domaine().nb_faces_elem();
 
  for (int elem=0; elem<nbelem; elem++)
    for (int i=0; i<nbsom_elem; i++)
      {
        // search for contributing faces, this is not optimal
        for (int f=0; f<nbface_elem; f++)
          {
            int face=elemfaces(elem,f);
            const double ratio = bidim_axi && xv_(face, 0) > 1e-10 ? xp_(elem, 0) / xv_(face, 0) : 1.0;
            for (int s=0; s<nbsom_face; s++)
              if (face_sommets_(face,s)==les_elems(elem,i))
 
                // search for contributing faces, this is not optimal
                for (int j=0; j<dimension; j++)
                  {
                    bij(elem,i,j)+=face_normales(face,j)*oriente_normale(face,elem) * ratio;
                  }
          }
      }
}
 
void Domaine_EF::calculer_Bij(DoubleTab& bij)
{
  int nbsom_face=nb_som_face();
  int nbelem=nb_elem_tot();
  int nbsom_elem=domaine().nb_som_elem();
  bij.resize(nbelem,nbsom_elem,dimension);
  const IntTab& les_elems=domaine().les_elems() ;
  if (nbsom_elem!=8)
    {
      calculer_Bij_gen(bij);
 
      bij/=nbsom_face;
      Bij_thilde_=bij;
      {
        //int nbelem=nb_elem_tot();
        //int nbsom_elem=domaine().nb_som_elem();
        const IntTab& elems=domaine().les_elems() ;
        //Cerr<<min(porosite_elem_)<<finl;abort();
        for (int elem=0; elem<nbelem; elem++)
          for (int i=0; i<nbsom_elem; i++)
            {
              int som_glob=elems(elem,i);;
              for (int j=0; j<dimension; j++)
                {
                  Bij_thilde_(elem,i,j)=bij(elem,i,j)*porosite_sommets_(som_glob);
                }
            }
 
        // essai
        if(0)
          {
            EFichier bthilde("Bthilde");
            for (int elem=0; elem<nbelem; elem++)
              for (int i=0; i<nbsom_elem; i++)
                for (int j=0; j<dimension; j++)
                  {
                    int ii,jj,ele;
                    //Cout<< elem<<" "<<i<<" "<<j<<" "<<Bij_thilde_(elem,i,j)<<" uuu ";
                    int i2=i;
                    if (i2==2) i2=3;
                    if (i2==3) i2=2;
                    if (i2==7) i2=6;
                    if (i2==6) i2=7;
                    Cout<< elem<<" "<<i<<" "<<j<<" "<<Bij_thilde_(elem,i2,j)<<" uuu ";
                    bthilde >> ele>> ii >> jj >> Bij_thilde_(elem,i2,j);
                    Cout<< Bij_thilde_(elem,i2,j)<<finl;
                    assert(ele==elem+1);
                    assert(ii==i+1);
                    assert(jj==j+1);
                  }
          }
      }
    }
  else
    {
      volumes_sommets_thilde_=0;
 
      Bij_thilde_.resize(nbelem,nbsom_elem,dimension);
      ArrOfInt filter(8) ;
 
 
      filter[0] = 0 ;
      filter[1] = 2 ;
      filter[2] = 3 ;
      filter[3] = 1 ;
      filter[4] = 4 ;
      filter[5] = 6 ;
      filter[6] = 7 ;
      filter[7] = 5 ;
      interface_INITGAUSS init_gauss;
      int npgau=27;
      // npgau=1;
      int nbnn=nbsom_elem;
      DoubleTab xgau(npgau,3),frgau(npgau,nbnn),dfrgau(npgau,nbnn,(int)3),poigau(npgau);
      DoubleVect volumes_sommets_(volumes_sommets_thilde_);
      int dim=3;
      //      init_gauss.compute(&dim,&nbnn,&npgau,xgau.addr(),frgau.addr(),dfrgau.addr(),poigau.addr());
      init_gauss.Compute(dim,nbnn,npgau,xgau,frgau,dfrgau,poigau);
      ArrOfInt num(nbnn);
      interface_CALCULBIJ CALCULBIJ;
      DoubleTab xl(dimension,nbnn),bijl(dimension,nbnn);
      DoubleTab detj(npgau),ajm1(npgau*9),aj(npgau*9),df(npgau*nbnn*3),iphi(nbnn);
      int ip=0;
      const DoubleTab& coord=domaine().coord_sommets();
      for (int elem=0; elem<nbelem; elem++)
        {
          for (int i=0; i<nbnn; i++)
            num[i]=les_elems(elem,filter[i]);
          for (int i=0; i<nbnn; i++)
            for (int d=0; d<dimension; d++)
              {
                xl(d,i)=coord(num[i],d);
              }
          num+=1;
          ip=0;
          double vol;
          int nbsom_tot=volumes_sommets_.size_totale();
          CALCULBIJ.Compute(nbnn,nbsom_tot, xl, num, bijl,
                            porosite_sommets_, ip,
                            npgau, xgau, frgau, dfrgau,
                            poigau,detj, ajm1, aj, df,vol,volumes_sommets_,iphi);
          if (!est_egal(volumes_(elem),vol,1e-5))
            {
              Cerr<<"Error volume "<<elem<< " vol: "<<volumes_(elem)<<" new "<<vol<<" delta "<<volumes_(elem)-vol;
              exit();
            }
          for (int i=0; i<nbnn; i++)
            {
              for (int d=0; d<dimension; d++)
                {
                  bij(elem,filter[i],d)=bijl(d,(i));
                  IPhi_(elem,filter[i])=iphi(i);
                }
            }
          ip=1;
 
          CALCULBIJ.Compute(nbnn, nbsom_tot, xl, num, bijl,
                            porosite_sommets_, ip,
                            npgau, xgau, frgau, dfrgau,
                            poigau,detj, ajm1, aj, df,volumes_thilde_(elem),volumes_sommets_thilde_,iphi);
 
          for (int i=0; i<nbnn; i++)
            {
              for (int d=0; d<dimension; d++)
                {
                  Bij_thilde_(elem,filter[i],d)=bijl(d,(i));
                  IPhi_thilde_(elem,filter[i])=iphi(i);
                }
 
            }
          //exit();
        }
    }
  if (0)
    {
      for (int elem=0; elem<nbelem; elem++)
        {
          for (int i=0; i<nbsom_elem; i++)
            {
              for (int d=0; d<dimension; d++)
                Cerr<<" BIJ "<<elem<<" "<< i<< " "<<d<<" "<<bij(elem,i,d)<<" "<<Bij_thilde_(elem,i,d)<<finl;
              Cerr<< " Iphi " << IPhi_(elem,i)<<" "<<IPhi_thilde_(elem,i)<<finl;
            }
          Cerr<<"vol elem "<<volumes_(elem)<< " "<<volumes_thilde_(elem)<<finl;
        }
      for (int elem=0; elem<volumes_sommets_thilde_.size_array(); elem++)
        Cerr<<elem<<" vol som "<< volumes_sommets_thilde_(elem)<<finl;
    }
  // verif
  Cerr<<"max/min/max_abs bij "<<mp_max_vect(bij)<< " "<<mp_min_vect(bij)<<" "<<mp_max_abs_vect(bij)<<finl;
  if (bidim_axi) return;
  int err=0;
  for (int elem=0; elem<nbelem; elem++)
    for (int j=0; j<dimension; j++)
      {
        double tot=0;
        for (int i=0; i<nbsom_elem; i++)
          tot+=bij(elem,i,j);
        if (!est_egal(tot,0))
          {
            err++;
            Cerr<<" pb elem "<< elem<<" direction "<<j<<" : "<<tot<<finl;
          }
        //Cerr<<"vol elem "<<volumes_(elem)<< " "<<volumes_thilde_(elem)<<finl;
      }
  if (err) exit();
 
}
void Domaine_EF::calculer_porosites_sommets()
{
 
 
  porosite_sommets_.resize(nb_som());
  porosite_sommets_=1;
}
void Domaine_EF::calculer_h_carre()
{
  // Computation of h_carre
  h_carre = 1.e30;
  h_carre_.resize(nb_faces());
  // Computation of face areas
  const DoubleVect& surfaces=face_surfaces();
  const int nb_faces_elem=domaine().nb_faces_elem();
  const int nbe=nb_elem();
  for (int num_elem=0; num_elem<nbe; num_elem++)
    {
      double surf_max = 0;
      for (int i=0; i<nb_faces_elem; i++)
        {
          double surf = surfaces(elem_faces(num_elem,i));
          surf_max = (surf > surf_max)? surf : surf_max;
        }
      double vol = volumes(num_elem)/surf_max;
      vol *= vol;
      h_carre_(num_elem) = vol;
      h_carre = ( vol < h_carre )? vol : h_carre;
    }
}
 
void Domaine_EF::calculer_volumes_entrelaces()
{
  //  Cerr << "les normales aux faces " << face_normales() << finl;
  // Compute the interleaved volumes;
  //  volumes_entrelaces_.resize(nb_faces());
  creer_tableau_faces(volumes_entrelaces_);
  int elem1,elem2;
  int nb_faces_elem=domaine().nb_faces_elem();
  int num_face;
  for (num_face=0; num_face<premiere_face_int(); num_face++)
    {
      elem1 = face_voisins_(num_face,0);
      volumes_entrelaces_[num_face] = volumes_[elem1]/nb_faces_elem;
    }
  for (num_face=premiere_face_int(); num_face<nb_faces(); num_face++)
    {
      elem1 = face_voisins_(num_face,0);
      elem2 = face_voisins_(num_face,1);
      volumes_entrelaces_[num_face] = (volumes_[elem1] + volumes_[elem2])
                                      /
                                      nb_faces_elem;
    }
  volumes_entrelaces_.echange_espace_virtuel();
}
 
void Domaine_EF::modifier_pour_Cl(const Conds_lim& conds_lim)
{
 
 
  //if (volumes_sommets_thilde_.size()!=nb_som())
  if (Bij_.nb_dim()==1)
    {
 
      calculer_volumes_sommets(conds_lim[0]->domaine_Cl_dis());
      Cerr << "le Domaine_EF a ete rempli avec succes" << finl;
 
      calculer_h_carre();
      // Compute porosities
 
      // vertex porosities are not used
      //calculer_porosites_sommets();
      calculer_Bij(Bij_);
 
 
 
    }
  Journal() << "Domaine_EF::Modifier_pour_Cl" << finl;
  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_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());
          int ndeb = 0;
          int nfin = la_front_dis.nb_faces_tot();
#ifndef NDEBUG
          int num_premiere_face = la_front_dis.num_premiere_face();
          int num_derniere_face = num_premiere_face+nfin;
#endif
          int nbr_faces_bord = la_front_dis.nb_faces();
          assert((nb_faces()==0)||(ndeb<nb_faces()));
          assert(nfin>=ndeb);
          int elem1,elem2,k;
          int face;
          // Modification of face_voisins_, face_normales_, volumes_entrelaces_ arrays.
          // Some normals are reoriented so that normals to periodic faces point
          // from face_voisins(face,0) toward face_voisins(face,1),
          // as is the case for internal faces.
 
          DoubleVect C1C2(dimension);
          double vol,psc=0;
 
          for (int ind_face=ndeb; ind_face<nfin; ind_face++)
            {
              //for ind_face
              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));
                  if (ind_face<nbr_faces_bord)
                    {
                      assert(faassociee>=num_premiere_face);
                      assert(faassociee<num_derniere_face);
                    }
 
                  elem1 = face_voisins_(face,0);
                  elem2 = face_voisins_(faassociee,0);
                  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;
                  psc = 0;
                  for (k=0; k<dimension; k++)
                    {
                      C1C2[k] = xv_(face,k) - xp_(face_voisins_(face,0),k);
                      psc += face_normales_(face,k)*C1C2[k];
                    }
 
                  if (psc < 0)
                    for (k=0; k<dimension; k++)
                      face_normales_(face,k) *= -1;
 
                  for (k=0; k<dimension; k++)
                    face_normales_(faassociee,k) = face_normales_(face,k);
                }
            }
        }
    }
 
  // PQ: 10/10/05: since periodic faces have a double contribution,
  //              the call to marquer_faces_double_contrib is made in this method
  //              in order to have access to conds_lim.
  Domaine_VF::marquer_faces_double_contrib(conds_lim);
}