Domaine_VF.h

/****************************************************************************
* 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.
*
*****************************************************************************/
 
#ifndef Domaine_VF_included
#define Domaine_VF_included
 
#include <Domaine_dis_base.h>
#include <Domaine_forward.h>
#include <TRUSTArrays.h>
#include <Front_VF.h>
 
#include <medcoupling++.h>
 
#ifdef MEDCOUPLING_
#include <MEDCouplingFieldTemplate.hxx>
#include <MEDCouplingCMesh.hxx>
using MEDCoupling::MEDCouplingCMesh;
using MEDCoupling::MCAuto;
using MEDCoupling::MEDCouplingFieldDouble;
#endif
 
class Domaine_Cl_dis_base;
 
/*! @brief class Domaine_VF
 *
 *  This abstract class contains the geometric information
 *  common to Finite Volume methods (VDF and VEF methods for example)
 *
 * @sa Domaine_dis_base
 */
class Domaine_VF : public Domaine_dis_base
{
  Declare_base(Domaine_VF);
public :
  virtual double face_normales(int face, int comp) const { return face_normales_(face,comp); }
  virtual DoubleTab& face_normales() { return face_normales_; }
  virtual const DoubleTab& face_normales() const { return face_normales_; }
  void calculer_face_surfaces(const DoubleVect& surfaces) { face_surfaces_ = surfaces; }
  virtual const DoubleVect& face_surfaces() const { return face_surfaces_; }
  virtual inline double face_surfaces(int i) const { return face_surfaces_(i); }
  virtual inline double surface(int i) const { return face_surfaces(i); }
 
  void discretiser() override;
  void discretiser_no_face() override;
  void typer_discretiser_ss_domaine(int i) override;
  void infobord();
  void info_elem_som();
  void marquer_faces_double_contrib(const Conds_lim&);
  virtual void typer_elem(Domaine&) {}
 
  virtual void remplir_face_voisins_fictifs(const Domaine_Cl_dis_base& ) ;
  virtual Faces* creer_faces();
  inline int nb_joints() const { return domaine().nb_joints(); }
  inline int premiere_face_int() const;
  inline int nb_faces() const;
  inline int nb_faces_tot() const;
  inline int nb_som_face() const;
  inline int nb_arete_face() const;
  inline int nb_faces_bord() const;
  inline int nb_faces_bord_tot() const;
  inline int premiere_face_bord() const;
  inline int nb_faces_internes() const;
 
  inline double xv(int num_face,int k) const { return xv_(num_face,k); }
  inline double xp(int num_elem,int k) const { return xp_(num_elem,k); }
  inline double xa(int num_arete,int k) const { return xa_(num_arete,k); }
 
  inline int face_numero_bord(int num_face) const;
  inline IntTab& face_numero_bord() { return face_numero_bord_; }
  inline const IntTab& face_numero_bord() const { return face_numero_bord_; }
  void remplir_face_numero_bord();
 
  inline ArrOfInt& est_face_bord() { return est_face_bord_; }
  inline const ArrOfInt& est_face_bord() const { return est_face_bord_; }
 
  inline virtual const IntVect& orientation() const;
  inline virtual int orientation(int ) const;
  inline virtual int orientation_si_definie(int) const;
 
  DoubleTab normalized_boundaries_outward_vector(int global_face_number, double scale_factor) const;
  inline DoubleTab& xv() { return xv_;}
  inline const DoubleTab& xv() const { return xv_;}
  inline DoubleTab& xp() { return xp_; }
  inline const DoubleTab& xp() const { return xp_; }
  inline DoubleTab& xa() { return xa_; }
  inline const DoubleTab& xa() const { return xa_; }
  inline DoubleVect& volumes_entrelaces() { return volumes_entrelaces_; }
  inline const DoubleVect& volumes_entrelaces() const { return volumes_entrelaces_; }
  inline double volumes_entrelaces(int num_face) const { return volumes_entrelaces_[num_face]; }
  inline const DoubleTab& volumes_entrelaces_dir() const { return volumes_entrelaces_dir_; }
  inline DoubleTab& volumes_entrelaces_dir() { return volumes_entrelaces_dir_; } // returns the array of interlaced volumes per side.
 
  inline const Joint& joint(int i) const { return domaine().joint(i); }
  inline Joint& joint(int i) { return domaine().joint(i); }
 
  inline Frontiere_dis_base& frontiere_dis(int ) override;
  inline const Frontiere_dis_base& frontiere_dis(int ) const override;
 
  inline int nb_frontiere_dis() const { return les_bords_.size(); }
  inline const Front_VF& front_VF(int i) const { return les_bords_[i]; } // returns the i-th discrete boundary.
  inline double volumes(int i) const { return volumes_[i]; }
  inline double inverse_volumes(int i) const { return inverse_volumes_[i]; }
  inline int face_voisins(int num_face,int i) const;
  inline int elem_faces(int i,int j) const;
  inline int face_sommets(int i,int j) const;
 
  inline DoubleVect& volumes() { return volumes_; }
  inline DoubleVect& inverse_volumes() { return inverse_volumes_; } // Array to optimize the code
  inline const DoubleVect& volumes() const { return volumes_; }
  inline const DoubleVect& inverse_volumes() const { return inverse_volumes_; } // Array to optimize the code
  inline IntTab& face_voisins() override;
  inline const IntTab& face_voisins() const override;
  inline const IntTab& face_voisins_fictifs() const { return face_voisins_fictifs_; }
  inline void face_voisins_reel_fictif(int face,int& el0,int& elf) const;
  inline IntTab& elem_faces();
  inline const IntTab& elem_faces() const;
  inline ArrOfInt& faces_doubles();
  inline const ArrOfInt& faces_doubles() const;
  inline IntTab& face_sommets() override;
  inline const IntTab& face_sommets() const override;
  inline IntTab& face_aretes() override;
  inline const IntTab& face_aretes() const override;
  void modifier_pour_Cl(const Conds_lim&) override;
 
  int numero_face_local(int face, int elem) const;
  inline int numero_sommet_local(int som, int elem) const;
 
  inline const IntTab& get_num_fac_loc() const { return num_fac_loc_; }
  inline int get_num_fac_loc(int, int) const;
  void construire_num_fac_loc();
 
  inline const ArrOfInt& ind_faces_virt_bord() const { return domaine().ind_faces_virt_bord(); }
  inline int est_une_face_virt_bord(int) const;
  inline int fbord(int f) const //returns the boundary face index of f if f is a boundary face, -1 otherwise
  {
    return f < premiere_face_int() ? f : f < nb_faces() ? -1 : ind_faces_virt_bord()[f - nb_faces()];
  }
 
  void construire_face_virt_pe_num();
  const IntTab& face_virt_pe_num() const;
 
  void creer_tableau_faces(Array_base&, RESIZE_OPTIONS opt = RESIZE_OPTIONS::COPY_INIT) const;
  void creer_tableau_aretes(Array_base&, RESIZE_OPTIONS opt = RESIZE_OPTIONS::COPY_INIT) const;
  void creer_tableau_faces_bord(Array_base&, RESIZE_OPTIONS opt = RESIZE_OPTIONS::COPY_INIT) const;
  const MD_Vector& md_vector_faces_bord() const { return md_vector_faces_front_; }
  const MD_Vector& md_vector_faces() const { return md_vector_faces_; }
  // Warning: if edges are not filled, md_vector_ is null
  const MD_Vector& md_vector_aretes() const { return md_vector_aretes_; }
 
  virtual const DoubleTab& xv_bord() const;
  DoubleTab calculer_xgr() const;
 
  virtual void get_position(DoubleTab& positions) const;
  virtual double compute_L1_norm(const DoubleVect& val_source, const bool basis_function, const int order) const;
  virtual double compute_L2_norm(const DoubleVect& val_source, const bool basis_function, const int order) const;
  virtual void compute_average(const DoubleVect& val_source, double& sum, double& volume, const bool basis_function, const int order) const;
  virtual void compute_average_porosity(const DoubleVect& val_source, const DoubleVect& porosity, double& sum, double& volume, const bool basis_function, const int order) const;
  virtual void get_nb_integ_points(IntTab& nelem) const;
  virtual void get_ind_integ_points(IntTab& nelem) const;
  virtual int get_max_nb_integ_points() const;
 
  //dot product (a - ma).(b - mb)
  inline double dot (const double *a, const double *b, const double *ma = nullptr, const double *mb = nullptr) const;
 
  //cross product
  inline std::array<double, 3> cross(int dima, int dimb, const double *a, const double *b, const double *ma = nullptr, const double *mb = nullptr) const;
 
  inline virtual double dist_norm(int ) const { Cerr << __func__ << " method should be overrided in a derived class !! " << finl; throw; }
  inline virtual double dist_norm_bord(int ) const { Cerr << __func__ << " method should be overrided in a derived class !! " << finl; throw; }
  inline virtual double dist_face_elem0(int ,int ) const { Cerr << __func__ << " method should be overrided in a derived class !! " << finl; throw; }
  inline virtual double dist_face_elem1(int ,int ) const { Cerr << __func__ << " method should be overrided in a derived class !! " << finl; throw; }
  inline virtual double dist_face_elem0_period(int ,int ,double ) const { Cerr << __func__ << " method should be overrided in a derived class !! " << finl; throw; }
  inline virtual double dist_face_elem1_period(int ,int ,double ) const { Cerr << __func__ << " method should be overrided in a derived class !! " << finl; throw; }
 
// Methods for computing and calling the distance to the nearest solid boundary; the input is the array of momentum equation BCs
  void init_dist_paroi_globale(const Conds_lim& conds_lim) override;
  const DoubleTab& normale_paroi_elem()  const {return n_y_elem_;}
  const DoubleTab& normale_paroi_faces() const {return n_y_faces_;}
 
  void build_mc_face_mesh() const;
  void build_mc_dual_mesh() const;
  inline int oriente_normale(int f, int e) const { return (face_voisins(f, 0) == e) ? 1 : -1; }
  inline const IntTab& get_face_voisins_dual() const { return face_dual_; }
  inline int get_face_voisins_dual(const int i, const int j) const { return face_dual_(i,j); }
 
#ifdef MEDCOUPLING_
  inline const MEDCouplingUMesh* get_mc_face_mesh() const;
  inline const MEDCouplingUMesh* get_mc_dual_mesh() const;
#endif
 
private:
  DoubleVect face_surfaces_;                // face surfaces
 
protected:
 
  DoubleVect volumes_;                          // element volumes
  DoubleVect inverse_volumes_;                  // inverse of element volumes
  DoubleVect volumes_entrelaces_;            // interlaced volumes for momentum integration
  DoubleTab volumes_entrelaces_dir_;        // interlaced volumes per side
  DoubleTab face_normales_;             // face normals
 
  inline double volume_entrelace_axi(double r_face, double r_elem, double axis_length) const;
 
  IntTab face_voisins_;                          // face/element connectivity
  IntTab face_voisins_fictifs_;           // face/fictitious element connectivity
  DoubleTab xp_;                            // element centres of gravity
  DoubleTab xv_;                            // face centres of gravity
  mutable DoubleTab xv_bord_;                       //xv_, but boundary faces only (created on demand)
 
  IntTab elem_faces_;                           // element/face connectivity
  IntTab face_sommets_;                           // face vertices
  IntTab face_aretes_;                           // face edges
  DoubleTab xa_;                            // edge centres of gravity
  IntTab face_numero_bord_;                     // face/boundary_number connectivity
 
  // Parallel descriptor for face arrays (size() == nb_faces())
  MD_Vector md_vector_faces_;
  // Same for boundary faces (size() == nb_faces_front())
  MD_Vector md_vector_faces_front_;
  // The one for edges
  MD_Vector md_vector_aretes_;
 
  VECT(Front_VF) les_bords_;
 
  IntTab num_fac_loc_;     // returns for a face its local number within the element
  //ArrOfInt faces_perio_;   // periodic faces (useful when looping from 0 to nb_faces_tot)
  ArrOfInt faces_doubles_; // faces with double contribution (periodic faces and shared items). Useful when looping from 0 to nb_faces for a subsequent reduction
  ArrOfInt est_face_bord_; // returns for a real or virtual face: 0 if internal, 1 if non-periodic boundary face, 2 if periodic boundary face
 
  // For each virtual face i with nb_faces_<=i<nb_faces_tot:
  // face_virt_pe_num_(i-nb_faces_,0) = number of the PE that owns the face
  // face_virt_pe_num_(i-nb_faces_,1) = local number of this face on the PE that owns it
  IntTab face_virt_pe_num_;
 
  DoubleTab n_y_elem_ ; // normal vector between the nearest boundary and the element
  DoubleTab n_y_faces_; // normal vector between the nearest boundary and the face
 
  int nb_elem_std_=-10;                     // number of standard elements
  int nb_faces_std_=-10;                    // number of standard faces
  IntVect rang_elem_non_std_;    // rang_elem_non_std_= -1 if the element is standard
  // rang_elem_non_std_= rank of the element in the arrays
  // relating to non-standard elements
 
  //
  // Dual mesh management:
  //
#ifdef MEDCOUPLING_
  ///! MEDCoupling version of the face domain - stored in Domaine_dis since faces are built here:
  mutable MCAuto<MEDCouplingUMesh> mc_face_mesh_;
  mutable MCAuto<MEDCouplingUMesh> mc_dual_mesh_;
  mutable bool mc_face_mesh_ready_ = false;
  mutable bool mc_dual_mesh_ready_ = false;
#endif
 
  mutable IntTab face_dual_; ///< For each face f, face_dual_(f, j) returns the element built on the left and right of the face in the dual mesh. Same sorting as face_voisins_
 
  void order_faces(Faces& les_faces);
  virtual void prepare_elem_non_std(Faces& les_faces);
  virtual void compute_sort_key(Faces& les_faces, IntTab& sort_key);
  void sort_along_zcurve(const Faces& les_faces, IntTab& sort_key) const;
  virtual void renumber_faces(Faces& les_faces, IntTab& sort_key);
 
  /*
   * XXX Elie Saikali
   *
   *  If requested, build a structured mesh of type MEDCouplingCMesh and the associated elem/face/node maps (CART -> TRUST) !
   *
   */
public:
  void build_map_mc_Cmesh(const bool with_faces) override;
 
 
#ifdef MEDCOUPLING_
private:
  MCAuto<MEDCouplingCMesh> mc_Cmesh_;
  std::vector<int> mc_Cmesh_elemCorrespondence_, mc_Cmesh_nodesCorrespondence_;
  std::vector<int> mc_Cmesh_facesXCorrespondence_, mc_Cmesh_facesYCorrespondence_, mc_Cmesh_facesZCorrespondence_;
  std::vector<double> mc_Cmesh_x_coords_, mc_Cmesh_y_coords_, mc_Cmesh_z_coords_;
 
  bool mc_Cmesh_ready_ = false, mc_Cmesh_with_faces_corr_ = false;
 
  void build_mc_Cmesh();
  void build_mc_Cmesh_nodesCorrespondence();
  void build_mc_Cmesh_correspondence(bool withFace);
 
  template <typename TYPE>
  TYPE Cmesh_error(const char * nom_funct) const
  {
    cerr << "Domaine_VF::" << nom_funct << " should not be called since the MEDCouplingCMesh of Domaine_VF is not yet filled !!!" << endl;
    Cerr << "Add the interpret Build_Map_to_Structured in your data file !!!" << finl;
    throw;
  }
 
public:
 
  inline const MEDCouplingCMesh* get_mc_CMesh() const
  {
    if (mc_Cmesh_ready_) return mc_Cmesh_;
    else return Cmesh_error<MEDCouplingCMesh*>(__func__);
  }
 
  const std::vector<int>& get_mc_Cmesh_elemCorrespondence() const
  {
    if (mc_Cmesh_ready_) return mc_Cmesh_elemCorrespondence_;
    else return Cmesh_error<std::vector<int>&>(__func__);
  }
 
  const std::vector<int>& get_mc_Cmesh_facesXCorrespondence() const
  {
    if (mc_Cmesh_ready_ && mc_Cmesh_with_faces_corr_) return mc_Cmesh_facesXCorrespondence_;
    else return Cmesh_error<std::vector<int>&>(__func__);
  }
 
  const std::vector<int>& get_mc_Cmesh_facesYCorrespondence() const
  {
    if (mc_Cmesh_ready_ && mc_Cmesh_with_faces_corr_) return mc_Cmesh_facesYCorrespondence_;
    else return Cmesh_error<std::vector<int>&>(__func__);
  }
 
  const std::vector<int>& get_mc_Cmesh_facesZCorrespondence() const
  {
    if (mc_Cmesh_ready_ && mc_Cmesh_with_faces_corr_) return mc_Cmesh_facesZCorrespondence_;
    else return Cmesh_error<std::vector<int>&>(__func__);
  }
 
  const std::vector<int>& get_mc_Cmesh_nodesCorrespondence() const
  {
    if (mc_Cmesh_ready_) return mc_Cmesh_nodesCorrespondence_;
    else return Cmesh_error<std::vector<int>&>(__func__);
  }
 
  const std::vector<double>& get_mc_Cmesh_x_coords() const
  {
    if (mc_Cmesh_ready_) return mc_Cmesh_x_coords_;
    else return Cmesh_error<std::vector<double>&>(__func__);
  }
 
  const std::vector<double>& get_mc_Cmesh_y_coords() const
  {
    if (mc_Cmesh_ready_) return mc_Cmesh_y_coords_;
    else return Cmesh_error<std::vector<double>&>(__func__);
  }
 
  const std::vector<double>& get_mc_Cmesh_z_coords() const
  {
    if (mc_Cmesh_ready_) return mc_Cmesh_z_coords_;
    else return Cmesh_error<std::vector<double>&>(__func__);
  }
 
  // Note: n_som i, not elem!
  int get_mc_Cmesh_ni() const
  {
    if (mc_Cmesh_ready_) return static_cast<int>(mc_Cmesh_x_coords_.size());
    else return Cmesh_error<int>(__func__);
  }
 
  // Note: n_som j, not elem!
  int get_mc_Cmesh_nj() const
  {
    if (mc_Cmesh_ready_) return static_cast<int>(mc_Cmesh_y_coords_.size());
    else return Cmesh_error<int>(__func__);
  }
 
  // Note: n_som k, not elem!
  int get_mc_Cmesh_nk() const
  {
    if (mc_Cmesh_ready_) return static_cast<int>(mc_Cmesh_z_coords_.size());
    else return Cmesh_error<int>(__func__);
  }
 
#endif
 
};
 
// Returns the local face number from a global face number and a local element index (0 or 1)
// Much faster method than Domaine_VF::numero_face_local(face,elem)
inline int Domaine_VF::get_num_fac_loc(int face,int voisin) const
{
  assert(voisin==0 || voisin==1);
  return num_fac_loc_(face,voisin);
}
 
inline int Domaine_VF::numero_sommet_local(int som, int elem) const
{
  int nse=domaine().nb_som_elem();
  const IntTab& les_elems=domaine().les_elems();
  for(int som_loc=0; som_loc<nse; som_loc++)
    {
      if(les_elems(elem, som_loc)==som)
        return som_loc;
    }
  return -1;
}
 
/*! @brief Returns the neighbouring element of num_face in direction i.
 *
 * i=0 : in the direction opposite to the axis orthogonal to face num_face.
 *  i=1 : in the direction of the axis orthogonal to face num_face.
 *  example:
 *
 *                 |   0     |     1    |
 *                       numface
 *  Returns -1 if the element does not exist (at the boundary).
 *
 */
inline int Domaine_VF::face_voisins(int num_face,int i) const
{
  return face_voisins_(num_face,i);
}
 
/*! @brief Returns the face-element connectivity array; see above.
 *
 */
inline IntTab& Domaine_VF::face_voisins()
{
  return face_voisins_;
}
 
/*! @brief See above.
 *
 */
inline const IntTab& Domaine_VF::face_voisins() const
{
  return face_voisins_;
}
 
/*! @brief Returns in el0 the index of the interior element; returns in elf the index of the fictitious element (-1 if it does not exist).
 *
 *
 *
 */
inline void  Domaine_VF::face_voisins_reel_fictif(int face,int& el0,int& elf) const
{
  assert(face<premiere_face_int());
  el0 = face_voisins(face,0);
  elf = face_voisins(face,1);
  if (el0==-1)
    {
      el0=elf;
      elf=face_voisins_fictifs_(face,0);
    }
  if (elf==-1) elf=face_voisins_fictifs_(face,1);
 
}
 
/*! @brief A face is internal if and only if it separates two elements.
 *
 * Returns the index of the first internal face.
 *
 */
inline int Domaine_VF::premiere_face_int() const
{
  return nb_faces_bord();
}
 
/*! @brief Returns the total number of faces.
 *
 */
inline int Domaine_VF::nb_faces() const
{
  return face_sommets_.dimension(0);
}
 
/*! @brief Returns the total number of faces.
 *
 * That is, real faces + virtual faces.
 *
 */
inline int Domaine_VF::nb_faces_tot() const
{
  return face_sommets_.dimension_tot(0);
}
 
 
/*! @brief Returns the number of vertices per face.
 *
 * Note: it is assumed that all faces have the same number of vertices,
 *  which excludes elements such as prisms.
 *
 */
inline int Domaine_VF::nb_som_face() const
{
  return face_sommets_.dimension(1);
}
 
/*! @brief Returns the number of edges per face.
 *
 *
 */
inline int Domaine_VF::nb_arete_face() const
{
  return face_aretes_.dimension(1);
}
 
/*! @brief Returns the number of faces on which boundary conditions are applied:
 *
 *  boundaries, connectors, plates.
 *
 */
inline int Domaine_VF::nb_faces_bord() const
{
  return domaine().nb_faces_frontiere();
}
 
/*! @brief Returns the index of the first face on which boundary conditions are applied:
 *
 *  boundaries, connectors, plates.
 *
 */
inline int Domaine_VF::premiere_face_bord() const
{
  return 0;
}
 
/*! @brief A face is internal if and only if it separates two elements.
 *
 * Returns the number of internal faces.
 *
 */
inline int Domaine_VF::nb_faces_internes() const
{
  return nb_faces()- nb_faces_bord();
}
 
/*! @brief Returns the index of the i-th face of element num_elem; the face numbering convention is
 *
 *  left to the responsibility of derived classes.
 *
 */
inline int Domaine_VF::elem_faces(int num_elem, int i) const
{
  return elem_faces_(num_elem, i);
}
 
/*! @brief Returns the element/face connectivity array.
 *
 */
inline IntTab& Domaine_VF::elem_faces()
{
  return elem_faces_;
}
 
/*! @brief See above.
 *
 */
inline const IntTab& Domaine_VF::elem_faces() const
{
  return elem_faces_;
}
 
/*! @brief Returns 1 for faces belonging to a periodic boundary or a shared item, 0 by default.
 *
 */
inline ArrOfInt& Domaine_VF::faces_doubles()
{
  return faces_doubles_;
}
 
/*! @brief See above.
 *
 */
inline const ArrOfInt& Domaine_VF::faces_doubles() const
{
  return faces_doubles_;
}
 
/*! @brief Returns the index of the i-th vertex of face num_face.
 *
 */
inline int Domaine_VF::face_sommets(int num_face, int i) const
{
  return face_sommets_(num_face, i);
}
 
/*! @brief Returns the face/vertex connectivity array.
 *
 */
inline IntTab& Domaine_VF::face_sommets()
{
  return face_sommets_;
}
 
/*! @brief See above.
 *
 */
inline const IntTab& Domaine_VF::face_sommets() const
{
  return face_sommets_;
}
 
/*! @brief Returns the face/edge connectivity array.
 *
 */
inline IntTab& Domaine_VF::face_aretes()
{
  return face_aretes_;
}
 
/*! @brief See above.
 *
 */
inline const IntTab& Domaine_VF::face_aretes() const
{
  return face_aretes_;
}
 
/*! @brief Returns the i-th discrete boundary.
 *
 */
inline const Frontiere_dis_base& Domaine_VF::frontiere_dis(int i) const
{
  return les_bords_[i];
}
 
/*! @brief Returns the i-th discrete boundary.
 *
 */
inline Frontiere_dis_base& Domaine_VF::frontiere_dis(int i)
{
  return les_bords_[i];
}
 
/*! @brief Returns the total number of faces (real + virtual) on which boundary conditions are applied:
 *
 *  boundaries, connectors, plates.
 *
 */
inline int Domaine_VF::nb_faces_bord_tot() const
{
  return md_vector_faces_bord()->get_nb_items_tot();
}
 
inline const IntVect& Domaine_VF::orientation() const
{
  Cerr<<"Try using orientation(num_face) instead of orientation() if the mesh is axis oriented."<< finl;
  exit();
  throw;
  //return orientation();
}
 
inline int Domaine_VF::orientation_si_definie(int num_face) const
{
  for (int dir=0; dir<dimension; dir++)
    {
      if (est_egal(std::fabs(face_normales(num_face, dir)), face_surfaces(num_face)))
        return dir;
    }
  return -1;
}
 
inline int Domaine_VF::orientation(int num_face) const
{
  for (int dir=0; dir<dimension; dir++)
    {
      if (est_egal(std::fabs(face_normales(num_face, dir)), face_surfaces(num_face)))
        return dir;
    }
  Cerr<<"Your mesh is not axis oriented."<<finl;
  Cerr<<"Try recoding using face_normales(). Example for g(face)=gravite[orientation[face]];" << finl;
  Cerr<<"Compute something like: g(face)=Sum(i=0;i<gravite.size();i++) of gravite[i]*face_normales(face,i)/surface(face);" << finl;
  exit();
  return -1;
}
 
/*! @brief Returns 1 if face is a virtual boundary face, 0 otherwise.
 *
 */
inline int Domaine_VF::est_une_face_virt_bord(int face) const
{
  if (face<nb_faces() || ind_faces_virt_bord()[face-nb_faces()]==-1)
    return 0;
  else
    return 1;
}
 
inline int Domaine_VF::face_numero_bord(int num_face) const
{
  assert(num_face < nb_faces());
  return face_numero_bord_(num_face);
}
 
/* dot product of two vectors */
inline double Domaine_VF::dot(const double *a, const double *b, const double *ma, const double *mb) const
{
  double res = 0;
  for (int i = 0; i < dimension; i++) res += (a[i] - (ma ? ma[i] : 0)) * (b[i] - (mb ? mb[i] : 0));
  return res;
}
 
inline double Domaine_VF::volume_entrelace_axi(double r_face, double r_elem, double axis_length) const
{
  const double dr = std::fabs(r_elem - r_face);
  const double r = (r_elem < r_face) ? r_elem : r_face;
  return 2. * M_PI * (r * dr + 0.5 * dr * dr) * axis_length;
}
 
/* cross product of two vectors (always 3D, even in 2D) */
inline std::array<double, 3> Domaine_VF::cross(int dima, int dimb, const double *a, const double *b, const double *ma, const double *mb) const
{
  std::array<double, 3> va = {{ 0, 0, 0 }}, vb = {{ 0, 0, 0 }}, res;
  for (int i = 0; i < dima; i++) va[i] = a[i] - (ma ? ma[i] : 0);
  for (int i = 0; i < dimb; i++) vb[i] = b[i] - (mb ? mb[i] : 0);
  for (int i = 0; i < 3; i++) res[i] = va[(i + 1) % 3] * vb[(i + 2) % 3] - va[(i + 2) % 3] * vb[(i + 1) % 3];
  // note : in 2D, the relevant results corresponding to the determinant is res[2]
  return res;
}
 
#ifdef MEDCOUPLING_
 
inline const MEDCouplingUMesh* Domaine_VF::get_mc_face_mesh() const
{
  if (!mc_face_mesh_ready_) build_mc_face_mesh();
  return mc_face_mesh_;
}
 
inline const MEDCouplingUMesh* Domaine_VF::get_mc_dual_mesh() const
{
  if (!mc_dual_mesh_ready_) build_mc_dual_mesh();
  return mc_dual_mesh_;
}
 
#endif
 
#endif /* Domaine_VF_included */