Domaine_PolyMAC_MPFA.cpp

/****************************************************************************
* Copyright (c) 2025, 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 <Echange_contact_PolyMAC_MPFA.h>
#include <Frottement_externe_impose.h>
#include <Linear_algebra_tools_impl.h>
#include <Frottement_global_impose.h>
#include <Champ_implementation_P1.h>
#include <Connectivite_som_elem.h>
#include <MD_Vector_composite.h>
#include <Dirichlet_homogene.h>
#include <MD_Vector_tools.h>
#include <Domaine_PolyMAC_MPFA.h>
#include <Domaine_Cl_PolyMAC_family.h>
#include <TRUSTTab_parts.h>
#include <Comm_Group_MPI.h>
#include <Quadrangle_VEF.h>
#include <communications.h>
#include <Hexaedre_VEF.h>
#include <Matrix_tools.h>
#include <unordered_map>
#include <Array_tools.h>
#include <TRUSTLists.h>
#include <Tetraedre.h>
#include <Rectangle.h>
#include <PE_Groups.h>
#include <Hexaedre.h>
#include <Triangle.h>
#include <EFichier.h>
#include <SFichier.h>
#include <Segment.h>
#include <Domaine.h>
#include <Scatter.h>
#include <EChaine.h>
#include <LireMED.h>
#include <Ecrire_MED.h>
#include <unistd.h>
#include <Lapack.h>
#include <numeric>
#include <vector>
#include <cfloat>
#include <tuple>
#include <cmath>
#include <cfenv>
#include <set>
#include <map>
 
Implemente_instanciable(Domaine_PolyMAC_MPFA, "Domaine_PolyMAC_MPFA", Domaine_PolyMAC_HFV);
 
Sortie& Domaine_PolyMAC_MPFA::printOn(Sortie& os) const { return Domaine_Poly_base::printOn(os); }
 
Entree& Domaine_PolyMAC_MPFA::readOn(Entree& is) { return Domaine_Poly_base::readOn(is); }
 
void Domaine_PolyMAC_MPFA::discretiser()
{
  /* skip the discretiser() of Domaine_PolyMAC_HFV to avoid initializing PolyMAC_HFV_V1 variables */
  Domaine_Poly_base::discretiser();
 
  //MD_vector for Champ_Face_PolyMAC_MPFA (faces + dimension * champ_p0)
  MD_Vector_composite mdc_ch_face;
  mdc_ch_face.add_part(md_vector_faces());
  mdc_ch_face.add_part(domaine().md_vector_elements(), dimension);
  mdv_ch_face.copy(mdc_ch_face);
}
 
/*! @brief Initializes the face stencils for gradient computation
 *
 * This method builds the connectivity stencils required for face-based gradient calculations
 * in the PolyMAC_MPFA discretization scheme. It establishes the relationship between faces
 * and all elements/boundary faces that are connected through shared vertices.
 *
 * The stencil for each face includes:
 * - All elements connected to the face through any of its vertices
 * - All boundary faces connected to the face through shared vertices
 * - The face's direct neighboring elements
 *
 * @note This method is called lazily - it only performs initialization if stencils
 *       haven't been built yet (checked via fsten_d.size())
 *
 * @note The method populates two main data structures:
 *       - fsten_d: Index array for stencil boundaries
 *       - fsten_eb: Flat array containing all stencil elements/boundary faces
 *
 * @details Algorithm overview:
 *          1. Build vertex-to-element connectivity (som_eb)
 *          2. Add boundary faces to vertex connectivity
 *          3. For each face, collect all elements/boundary faces connected via vertices
 *          4. Store results in compressed sparse format
 *
 * @post After execution:
 *       - fsten_d(f) gives the starting index in fsten_eb for face f's stencil
 *       - fsten_d(f+1) gives the ending index (exclusive) for face f's stencil
 *       - fsten_eb contains the actual element/boundary face indices
 *       - Boundary faces are represented as (ne_tot + face_index)
 *
 * @note Memory is optimized using CRIMP() to reduce storage overhead
 *
 * @note Only processes faces that have at least one valid neighbor:
 *       - Internal faces: both neighbors must be valid
 *       - Boundary faces: at least one neighbor must be valid
 *
 * @see fgrad() which uses these stencils for gradient computation
 */
void Domaine_PolyMAC_MPFA::init_stencils() const
{
  if (fsten_d.size())
    return;
 
  const IntTab& f_s = face_sommets(), &f_e = face_voisins(), &e_s = domaine().les_elems();
  const int ne_tot = nb_elem_tot(), ns_tot = domaine().nb_som_tot();
  fsten_d.resize(1);
 
  /* vertex -> element / boundary face connectivity */
  std::vector<std::set<int>> som_eb(ns_tot);
  for (int e = 0; e < nb_elem_tot(); e++)
    for (int i = 0; i < e_s.dimension(1); i++)
      {
        const int s = e_s(e, i);
        if (s < 0) continue;
 
        som_eb[s].insert(e);
      }
 
  for (int f = 0; f < nb_faces_tot(); f++)
    if (fbord(f) >= 0)
      for (int i = 0; i < f_s.dimension(1); i++)
        {
          const int s = f_s(f, i);
          if (s < 0) continue;
 
          som_eb[s].insert(ne_tot + f);
        }
 
  std::set<int> f_eb; //vertices of face f, elements connected to e via vertices, vertices / faces connected through a common face
  for (int f = 0; f < nb_faces_tot(); fsten_d.append_line(fsten_eb.size()), f++)
    if (f_e(f, 0) >= 0 && (fbord(f) >= 0 || f_e(f, 1) >= 0))
      {
        /* connectivity through a vertex of f */
        f_eb.clear();
        for (int i = 0; i < f_s.dimension(1); i++)
          {
            const int s = f_s(f, i);
            if (s < 0) continue;
 
            for (auto &&el : som_eb[s])
              f_eb.insert(el);
          }
 
        /* fill */
        for (auto &&el : f_eb)
          fsten_eb.append_line(el);
      }
 
  CRIMP(fsten_d);
  CRIMP(fsten_eb);
}
 
/*! @brief Computes field gradient interpolation at faces while preserving constant fields
 *
 * This method computes the interpolation [n_f.grad T]_f (if nu_grad = 0) or [n_f.nu.grad T]_f
 * while exactly preserving fields satisfying [nu grad T]_e = constant.
 * It uses MPFA (Multi-Point Flux Approximation) methods with three strategies:
 *       - Attempt 0: Standard MPFA-O
 *       - Attempt 1: MPFA-eta: variant of MPFA-O where auxiliary variables are not in the barycenter of the face
 *       - Attempt 2: Symmetric MPFA if previous attempts fail (coercive but not always consistent, especially on complex meshes)
 *
 * @param N             Number of field components
 * @param is_p          Pressure field indicator (1 for pressure, 0 otherwise)
 *                      Controls Neumann/Dirichlet inversion for boundary conditions
 * @param cls           Domain boundary conditions
 * @param fcl           Boundary condition data
 *                      fcl(f, 0/1/2): (BC type, BC index, index within BC)
 *                      See Champ_{P0,Face}_PolyMAC_MPFA for format
 * @param nu            Element diffusivity (optional, can be nullptr)
 *                      Array nu(e, n, ..) for element e and component n
 * @param som_ext       List of vertices to exclude from processing (optional)
 *                      Example: direct treatment of Echange_Contact in Op_Diff_PolyMAC_MPFA_Elem
 * @param virt          Virtual faces indicator (1 to include, 0 otherwise)
 * @param full_stencil  Complete stencil indicator (1 for full dimensioning)
 *
 * @param[out] phif_d   Start/end indices in phif_{e,c} / phif_{pe,pc}
 *                      phif_d(f, 0/1): flux indices at face f in interval
 *                      [phif_d(f, 0/1), phif_d(f + 1, 0/1)[
 * @param[out] phif_e   Element indices in stencil for each contribution
 *                      phif_e(i): element index for i-th contribution
 * @param[out] phif_c   Stencil coefficients
 *                      phif_c(i, n, c): coefficient for element i, component n, contribution c
 *                      Contains local indices/coefficients (without Echange_contact)
 *                      and diagonal terms (independent components)
 *
 * @note The method checks positivity of bilinear form eigenvalues to ensure scheme stability and choose the MPFA method accordingly
 *
 * @note The method returns also a percentage of which MPFA method is used. Be careful of the validity of the solution if the percentage of MPFA-sym is high
 *
 * @note Special handling for different boundary condition types:
 *       - Dirichlet/Neumann
 *       - Imposed global/external friction
 *       - Imposed global/external exchange
 *
 */
void Domaine_PolyMAC_MPFA::fgrad(int N, int is_p, const Conds_lim& cls, const IntTab& fcl, const DoubleTab *nu, const IntTab *som_ext,
                                 int virt, int full_stencil, IntTab& phif_d, IntTab& phif_e, DoubleTab& phif_c) const
{
#ifdef _COMPILE_AVEC_PGCC_AVANT_22_7
  Cerr << "Internal error with nvc++: Internal error: read_memory_region: not all expected entries were read." << finl;
  Process::exit();
#else
  const IntTab& f_e = face_voisins(), &e_f = elem_faces(), &f_s = face_sommets();
  const DoubleTab& nf = face_normales(), &xs = domaine().coord_sommets();
  const DoubleVect& fs = face_surfaces(), &vf = volumes_entrelaces();
  const Static_Int_Lists& s_e = som_elem();
  int i, i_s, j, k, l, e, f, s, sb, n_f, n_m, n_ef, n_e, n_eb, m, n, ne_tot = nb_elem_tot(), sgn, nw, infoo=-1, d, db,
                                                                     D = dimension, rk=-1, nl, nc, un = 1, il, ok, essai;
 
  unsigned long ll;
 
  double x, eps_g = 1e-6, eps = 1e-10, i3[3][3] = { { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 }}, fac[3], vol_s;
 
  init_stencils();
  phif_e.resize(0);
  phif_c.resize(fsten_eb.dimension(0), N);
  phif_c = 0;
 
  std::vector<int> s_eb, s_f; //lists of elements/boundary items, faces around the vertex
  std::vector<double> surf_fs, vol_es; //partial surfaces of faces connected to the vertex (same order as s_f)
  std::vector<std::array<std::array<double, 3>,2>> vec_fs;//for each facet, a basis of (D-1) vectors spanning it
  std::vector<std::vector<int>> se_f; /* se_f[i][.]: faces connected to the i-th element connected to vertex s */
  DoubleTrav M, B, X, Ff, Feb, Mf, Meb, W(1), x_fs, A, S; //system M.(grad u) = B in each element, flux at face Ff.u_fs + Feb.u_eb, equations Mf.u_fs = Meb.u_eb
  IntTrav piv, ctr[3];
  for (i = 0; first_fgrad_ && i < 3; i++) domaine().creer_tableau_sommets(ctr[i]);
 
  /* contributions at vertices: skipping those in som_ext */
  for (i_s = 0; i_s <= (som_ext ? som_ext->size() : 0); i_s++)
    for (s = (som_ext && i_s ? (*som_ext)(i_s - 1) + 1 : 0); s < (som_ext && i_s < som_ext->size() ? (*som_ext)(i_s) : (virt ? nb_som_tot() : nb_som())); s++)
      {
        /* elements connected to s: from som_elem (already sorted) */
        for (s_eb.clear(), n_e = 0; n_e < s_e.get_list_size(s); n_e++) s_eb.push_back(s_e(s, n_e));
        /* faces and their partial surfaces */
        for (s_f.clear(), surf_fs.clear(), vec_fs.clear(), se_f.resize(std::max(int(se_f.size()), n_e)), i = 0, ok = 1; i < n_e; i++)
          for (se_f[i].clear(), e = s_eb[i], j = 0; j < e_f.dimension(1) && (f = e_f(e, j)) >= 0; j++)
            {
              for (k = 0, sb = 0; k < f_s.dimension(1) && (sb = f_s(f, k)) >= 0; k++)
                if (sb == s) break;
              if (sb != s) continue; /* face of e not connected to s -> skip */
              if (fbord(f) >= 0) s_eb.insert(std::lower_bound(s_eb.begin(), s_eb.end(), ne_tot + f), ne_tot + f); //if f is a boundary face, add its index to s_eb
              else ok &= (f_e(f, 0) >= 0 && f_e(f, 1) >= 0); //if f is internal, upstream/downstream must be present
              se_f[i].push_back(f); //faces connected to both e and s
              if ((ll = std::lower_bound(s_f.begin(), s_f.end(), f) - s_f.begin()) == s_f.size() || s_f[ll] != f) /* if f is not in s_f, add it */
                {
                  s_f.insert(s_f.begin() + ll, f); //add face to s_f
                  if (D < 3) surf_fs.insert(surf_fs.begin() + ll, fs(f) / 2), vec_fs.insert(vec_fs.begin() + ll, {{{ xs(s, 0) - xv_(f, 0), xs(s, 1) - xv_(f, 1), 0}, { 0, 0, 0 }}}); //2D -> straightforward
                  else for (surf_fs.insert(surf_fs.begin() + ll, 0), vec_fs.insert(vec_fs.begin() + ll, {{{ 0, 0, 0}, {0, 0, 0 }}}), m = 0; m < 2; m++) //3D -> two sub-triangles
                  {
                    if (m == 1 || k > 0) sb = f_s(f, m ? (k + 1 < f_s.dimension(1) && f_s(f, k + 1) >= 0 ? k + 1 : 0) : k - 1); //next vertex (m = 1) or previous with k > 0 -> straightforward
                    else for (n = f_s.dimension(1) - 1; (sb = f_s(f, n)) == -1; ) n--; //previous vertex with k = 0 -> search from the end
                    auto v = cross(D, D, &xs(s, 0), &xs(sb, 0), &xv_(f, 0), &xv_(f, 0));//cross product (xs - xf)x(xsb - xf)
                    if (fs(f) > 0) surf_fs[ll] += std::fabs(dot(&v[0], &nf(f, 0))) / fs(f) / 4; //surface to add
                    for (d = 0; d < D; d++) vec_fs[ll][m][d] = (xs(s, d) + xs(sb, d)) / 2 - xv_(f, d); //vector face -> edge
                  }
                }
            }
        if (!ok) continue; //at least one missing neighbor
        n_eb = (int)s_eb.size(), n_f = (int)s_f.size();
 
        /* convert se_f to indices in s_f */
        for (i = 0; i < n_e; i++)
          for (j = 0; j < (int) se_f[i].size(); j++) se_f[i][j] = (int)(std::lower_bound(s_f.begin(), s_f.end(), se_f[i][j]) - s_f.begin());
        for (vol_es.resize(n_e), vol_s = 0, i = 0; i < n_e; vol_s += vol_es[i], i++)
          for (e = s_eb[i], vol_es[i] = 0, j = 0; j < (int) se_f[i].size(); j++)
            {
              f = s_f[k = se_f[i][j]];
              if (fs(f) > 0) vol_es[i] += surf_fs[k] * std::fabs(dot(&xp_(e, 0), &nf(f, 0), &xv_(f, 0))) / fs(f) / D;
            }
 
        for (essai = 0; essai < 3; essai++) /* attempt 0: MPFA-O -> attempt 1: MPFA-O with mobile x_fs -> attempt 2: symmetric MPFA (coercive but not very consistent) */
          {
            if (essai == 1) /* attempt 1: choose x_fs continuity points at facets to obtain a symmetric scheme */
              {
                /* linear system */
                for (M.resize(N, nc = (D - 1) * n_f, nl = D * (D - 1) / 2 * n_e), B.resize(N, n_m = std::max(nc, nl)), M = 0, B = 0, i = 0, il = 0; i < n_e; i++)
                  for (d = 0; d < D; d++)
                    for (db = 0; db < d; db++, il++)
                      for (e = s_eb[i], j = 0; j < (int) se_f[i].size(); j++)
                        for (sgn = e == f_e(f = s_f[k = se_f[i][j]], 0) ? 1 : -1, n = 0; n < N; n++)
                          {
                            const double inv_fs = fs(f) > 0 ? 1. / fs(f) : 0.;
                            for (l = 0; l < D; l++) fac[l] = sgn * nu_dot(nu, e, n, &nf(f, 0), i3[l]) * surf_fs[k] * inv_fs / vol_es[i]; //outward lambda_e nf vector * common factor
                            B(n, il) += fac[d] * (xv_(f, db) - xp_(e, db)) - fac[db] * (xv_(f, d) - xp_(e, d)); //right-hand side
                            for (l = 0; l < D - 1; l++) M(n, (D - 1) * k + l, il) += fac[db] * vec_fs[k][l][d] - fac[d] * vec_fs[k][l][db]; //matrix
                          }
 
                /* solve -> DGELSY */
                nw = -1, piv.resize(nc), F77NAME(dgelsy)(&nl, &nc, &un, &M(0, 0, 0), &nl, &B(0, 0), &n_m, &piv(0), &eps_g, &rk, &W(0), &nw, &infoo);
                for (W.resize(nw = (int)std::lrint(W(0))), n = 0; n < N; n++) piv = 0, F77NAME(dgelsy)(&nl, &nc, &un, &M(n, 0, 0), &nl, &B(n, 0), &n_m, &piv(0), &eps_g, &rk, &W(0), &nw, &infoo);
                /* x_fs = xf + corrections (clamped to [0, 0.5]) */
                for (x_fs.resize(N, n_f, D), n = 0; n < N; n++)
                  for (i = 0; i < n_f; i++)
                    for (f = s_f[i], d = 0; d < D; d++)
                      for (x_fs(n, i, d) = xv_(f, d), k = 0; k < D - 1; k++) x_fs(n, i, d) += std::min(std::max(B(n, (D - 1) * i + k), 0.), 0.5) * vec_fs[i][k][d];
              }
 
            /* per-cell gradients as functions of (u_eb, u_fs), flux F = Ff.u_fs + Feb.u_eb, and system Mf.u_fs = Feb.u_eb */
            Ff.resize(n_f, n_f, N), Feb.resize(n_f, n_eb, N), Mf.resize(N, n_f, n_f), Meb.resize(N, n_eb, n_f);
            for (Ff = 0, Feb = 0, Mf = 0, Meb = 0, i = 0; i < n_e; i++)
              for (e = s_eb[i], M.resize(n_ef = (int)se_f[i].size(), D), B.resize(D, n_m = std::max(D, n_ef)), X.resize(n_ef, D), piv.resize(n_ef), n = 0; n < N; n++)
                {
                  if (essai < 2) /* attempts 0 and 1: consistent gradient given by (u_e, (u_fs)_{f adj e, s})*/
                    {
                      /* gradient in (e, s) -> matrix / right-hand side M.x = B of the system (grad u)_i = sum_f b_{fi} (x_fs_i - x_e), where x_fs is the continuity point of u_fs */
                      for (j = 0; j < n_ef; j++)
                        for (f = s_f[k = se_f[i][j]], d = 0; d < D; d++) M(j, d) = (essai ? x_fs(n, k, d) : xv_(f, d)) - xp_(e, d);
                      for (B = 0, d = 0; d < D; d++) B(d, d) = 1;
                      nw = -1, piv = 0, F77NAME(dgelsy)(&D, &n_ef, &D, &M(0, 0), &D, &B(0, 0), &n_m, &piv(0), &eps_g, &rk, &W(0), &nw, &infoo);
                      W.resize(nw = (int)std::lrint(W(0))), F77NAME(dgelsy)(&D, &n_ef, &D, &M(0, 0), &D, &B(0, 0), &n_m, &piv(0), &eps_g, &rk, &W(0), &nw, &infoo);
                      for (j = 0; j < n_ef; j++)
                        for (d = 0; d < D; d++) X(j, d) = B(d, j); /* transposed to enable nu_dot */
                    }
                  else for (j = 0; j < n_ef; j++)
                      for (sgn = e == f_e(f = s_f[k = se_f[i][j]], 0) ? 1 : -1, d = 0; d < D; d++) /* attempt 2: inconsistent gradient */
                        {
                          const double inv_fs = fs(f) > 0 ? 1. / fs(f) : 0.;
                          X(j, d) = surf_fs[k] / vol_es[i] * sgn * nf(f, d) * inv_fs;
                        }
 
                  /* flux and equation. Note: complex BCs for scalar equations are handled directly in Op_Diff_PolyMAC_MPFA_Elem */
                  for (j = 0; j < n_ef; j++)
                    {
                      k = se_f[i][j], f = s_f[k], sgn = e == f_e(f, 0) ? 1 : -1; //face and its index
                      const Cond_lim_base *cl = fcl(f, 0) ? &cls[fcl(f, 1)].valeur() : nullptr; //pointer to BC if on a boundary
                      int is_dir = cl && (is_p ? sub_type(Neumann, *cl) : sub_type(Dirichlet, *cl) || sub_type(Dirichlet_homogene, *cl)); //is it Dirichlet?
                      for (l = 0; l < n_ef; l++)
                        {
                          const double inv_fs = fs(f) > 0 ? 1. / fs(f) : 0.;
                          x = sgn * nu_dot(nu, e, n, &nf(f, 0), &X(l, 0)) * surf_fs[k] * inv_fs; //contribution to the flux
                          if (sgn > 0) Ff(k, se_f[i][l], n) += x, Feb(k, i, n) -= x; //flux upstream->downstream
                          if (!is_dir) Mf(n, se_f[i][l], k) += x, Meb(n, i, k) += x; //equation on u_fs (except if Dirichlet BC)
                        }
                      if (!cl) continue; //nothing on the other side
                      else if (is_dir) Mf(n, k, k) = Meb(n, (int)(std::find(s_eb.begin(), s_eb.end(), ne_tot + f) - s_eb.begin()), k) = 1; //Dirichlet -> equation u_fs = u_b
                      else if (is_p ? !is_dir : sub_type(Neumann, *cl)) //Neumann -> add boundary flux
                        Meb(n, (int)(std::find(s_eb.begin(), s_eb.end(), ne_tot + f) - s_eb.begin()), k) += surf_fs[k];
                      else if (sub_type(Frottement_global_impose, *cl)) //Frottement_global_impose -> flux =  - coeff * v_e
                        Meb(n, i, k) -= surf_fs[k] * ((nu) ? ref_cast(Frottement_global_impose, *cl).coefficient_frottement(fcl(f, 2), n) : ref_cast(Frottement_global_impose, *cl).coefficient_frottement_grad(fcl(f, 2), n) ) ;
                      else if (sub_type(Frottement_externe_impose, *cl)) //Frottement_externe_impose -> flux =  - coeff * v_f
                        Mf(n, k, k) += surf_fs[k] * ((nu) ? ref_cast(Frottement_externe_impose, *cl).coefficient_frottement(fcl(f, 2), n) : ref_cast(Frottement_externe_impose, *cl).coefficient_frottement_grad(fcl(f, 2), n) );
                      else if (sub_type(Echange_impose_base, *cl)) //Echange_impose_base -> flux =  - h * (T_{e,f} - T_ext)
                        {
                          double h = (nu) ? ref_cast(Echange_impose_base, *cl).h_imp(fcl(f, 2), n) : ref_cast(Echange_impose_base, *cl).h_imp_grad(fcl(f, 2), n) ;
                          Meb(n, (int)(std::find(s_eb.begin(), s_eb.end(), ne_tot + f) - s_eb.begin()), k) += surf_fs[k] * h; //h * T_ext contribution
                          if (sub_type(Echange_externe_impose, *cl)) Mf(n, k, k) += surf_fs[k] * h; //Echange_externe_impose: h * T_f contribution
                          else Meb(n, i, k) -= surf_fs[k] * h; //Echange_global_impose: h * T_e contribution
                        }
                    }
                }
            /* solve Mf.u_fs = Meb.u_eb: DGELSY, just in case */
            nw = -1, piv.resize(n_f), F77NAME(dgelsy)(&n_f, &n_f, &n_eb, &Mf(0, 0, 0), &n_f, &Meb(0, 0, 0), &n_f, &piv(0), &eps, &rk, &W(0), &nw, &infoo);
            for (W.resize(nw = (int)std::lrint(W(0))), n = 0; n < N; n++)
              piv = 0, F77NAME(dgelsy)(&n_f, &n_f, &n_eb, &Mf(n, 0, 0), &n_f, &Meb(n, 0, 0), &n_f, &piv(0), &eps, &rk, &W(0), &nw, &infoo);
 
            /* substitute into Feb */
            for (i = 0; i < n_f; i++)
              for (j = 0; j < n_eb; j++)
                for (n = 0; n < N; n++)
                  for (k = 0; k < n_f; k++)
                    Feb(i, j, n) += Ff(i, k, n) * Meb(n, j, k);
 
            /* A: bilinear form */
            if (essai == 2) break;//not needed for VFSYM
            for (A.resize(N, n_e, n_e), A = 0, i = 0; i < n_e; i++)
              for (e = s_eb[i], j = 0; j < (int) se_f[i].size(); j++)
                for (sgn = e == f_e(f = s_f[k = se_f[i][j]], 0) ? 1 : -1, l = 0; l < n_e; l++)
                  for (n = 0; n < N; n++)
                    A(n, i, l) -= sgn * Feb(k, l, n);
            /* symmetrization */
            for (n = 0; n < N; n++)
              for (i = 0; i < n_e; i++)
                for (j = 0; j <= i; j++) A(n, i, j) = A(n, j, i) = (A(n, i, j) + A(n, j, i)) / 2;
            /* smallest eigenvalue: DSYEV */
            nw = -1, F77NAME(DSYEV)("N", "U", &n_e, &A(0, 0, 0), &n_e, S.addr(), &W(0), &nw, &infoo);
            for (W.resize(nw = (int)std::lrint(W(0))), S.resize(n_e), n = 0, ok = 1; n < N; n++)
              F77NAME(DSYEV)("N", "U", &n_e, &A(n, 0, 0), &n_e, &S(0), &W(0), &nw, &infoo), ok &= S(0) > -1e-8 * vol_s;
            if (ok) break; //so that "essai" has the correct value on exit
          }
        if (first_fgrad_) ctr[essai](s) = 1;
 
        /* store into phif_c */
        for (i = 0; i < n_f; i++)
          for (f = s_f[i], j = 0; j < n_eb; j++)
            for (k = (int)(std::lower_bound(fsten_eb.addr() + fsten_d(f), fsten_eb.addr() + fsten_d(f + 1), s_eb[j]) - fsten_eb.addr()), n = 0; n < N; n++)
              if (fs(f) > 0) phif_c(k, n) += Feb(i, j, n) / fs(f);
      }
 
 
  /* stencil simplification */
 
  int skip;
  DoubleTrav scale(N);
  for (phif_d.resize(1), phif_d = 0, phif_e.resize(0), f = 0, i = 0; f < nb_faces_tot(); f++, phif_d.append_line(i))
    if (fbord(f) >= 0 || (f_e(f, 0) >= 0 && f_e(f, 1) >= 0))
      {
        for (n = 0; n < N; n++)
          if (fs(f) > 0) scale(n) = nu_dot(nu, f_e(f, 0), n, &nf(f, 0), &nf(f, 0)) / (fs(f) * vf(f)); //order of magnitude of the coefficients
        for (j = fsten_d(f); j < fsten_d(f + 1); j++)
          {
            for (skip = !full_stencil && fsten_eb(j) != f_e(f, 0), n = 0; n < N; n++) skip &= std::fabs(phif_c(j, n)) < 1e-8 * scale(n); //what threshold to use here?
            if (skip) continue;
            for (n = 0; n < N; n++) phif_c(i, n) = phif_c(j, n);
            phif_e.append_line(fsten_eb(j)), i++;
          }
      }
  /* counting */
  if (!first_fgrad_) return;
  double count[3] = { mp_somme_vect_as_double(ctr[0]), mp_somme_vect_as_double(ctr[1]), mp_somme_vect_as_double(ctr[2]) }, tot = count[0] + count[1] + count[2];
  if (tot > 1.0e-4)
    Cerr << domaine().le_nom() << "::fgrad(): " << 100. * count[0] / tot << "% MPFA-O "
         << 100. * count[1] / tot << "% MPFA-O(h) " << 100. * count[2] / tot << "% MPFA-SYM" << finl;
  first_fgrad_ = 0;
#endif
}