Champ_Face_PolyMAC_CDO.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 <Linear_algebra_tools_impl.h>
#include <Connectivite_som_elem.h>
#include <Dirichlet_homogene.h>
#include <Domaine_Cl_PolyMAC_family.h>
#include <Champ_Fonc_reprise.h>
#include <Champ_Face_PolyMAC_CDO.h>
#include <Schema_Temps_base.h>
#include <Champ_Uniforme.h>
#include <Domaine_PolyMAC_CDO.h>
 
#include <TRUSTTab_parts.h>
#include <Probleme_base.h>
#include <Equation_base.h>
#include <Matrix_tools.h>
#include <Domaine_VF.h>
#include <Dirichlet.h>
#include <Symetrie.h>
#include <EChaine.h>
#include <Domaine.h>
#include <array>
#include <cmath>
 
Implemente_instanciable(Champ_Face_PolyMAC_CDO, "Champ_Face_PolyMAC_CDO|Champ_Face_PolyMAC", Champ_Face_base);
 
Sortie& Champ_Face_PolyMAC_CDO::printOn(Sortie& os) const { return os << que_suis_je() << " " << le_nom(); }
 
Entree& Champ_Face_PolyMAC_CDO::readOn(Entree& is) { return is; }
 
int Champ_Face_PolyMAC_CDO::fixer_nb_valeurs_nodales(int n)
{
  // use the same kind of code as in Champ_Fonc_P0_base, but retrieve the number of faces instead of the number of elements
  // note: the code used in Champ_Fonc_P0_base::fixer_nb_valeurs_nodales() to retrieve the discrete domain is somewhat surprising...
 
  const Champ_Inc_base& self = ref_cast(Champ_Inc_base, *this);
  const Domaine_PolyMAC_CDO& domaine = ref_cast(Domaine_PolyMAC_CDO,self.domaine_dis_base());
 
  assert(n == domaine.nb_faces() + (dimension < 3 ? domaine.nb_som() : domaine.domaine().nb_aretes()));
 
  // Problem: nb_comp equals 2 but we only want one dimension !!!
  // HACK:
  int old_nb_compo = nb_compo_;
  nb_compo_ = 1;
  creer_tableau_distribue(domaine.mdv_faces_aretes);
  nb_compo_ = old_nb_compo;
  return n;
}
 
Champ_base& Champ_Face_PolyMAC_CDO::affecter_(const Champ_base& ch)
{
  const DoubleTab& v = ch.valeurs();
  DoubleTab_parts parts(valeurs());
  DoubleTab& val = parts[0]; //velocity part
  const Domaine_VF& domaine_PolyMAC_CDO = ref_cast( Domaine_VF,le_dom_VF.valeur());
  int nb_faces = domaine_PolyMAC_CDO.nb_faces();
  const DoubleVect& surface = domaine_PolyMAC_CDO.face_surfaces();
  const DoubleTab& normales = domaine_PolyMAC_CDO.face_normales();
 
  if (sub_type(Champ_Uniforme,ch))
    {
      for (int num_face=0; num_face<nb_faces; num_face++)
        {
          double vn=0;
          for (int dir=0; dir<dimension; dir++)
            vn+=v(0,dir)*normales(num_face,dir);
 
          vn/=surface(num_face);
          val(num_face) = vn;
        }
    }
  else if (sub_type(Champ_Fonc_reprise, ch))
    for (int f = 0; f < nb_faces; f++)
      val(f) = ch.valeurs()(f);
  else
    {
      //      int ndeb_int = domaine_PolyMAC_CDO.premiere_face_int();
      //      const IntTab& face_voisins = domaine_PolyMAC_CDO.face_voisins();
      const DoubleTab& xv=domaine_PolyMAC_CDO.xv();
      DoubleTab eval(val.dimension_tot(0),dimension);
      ch.valeur_aux(xv,eval);
      for (int num_face=0; num_face<nb_faces; num_face++)
        {
          double vn=0;
          for (int dir=0; dir<dimension; dir++)
            vn+=eval(num_face,dir)*normales(num_face,dir);
 
          vn/=surface(num_face);
          val(num_face) = vn;
        }
    }
  return *this;
}
 
DoubleVect& Champ_Face_PolyMAC_CDO::valeur_a_elem(const DoubleVect& position, DoubleVect& result, int poly) const
{
  Cerr << "Champ_Face_PolyMAC_CDO::" <<__func__ << " is not coded !" << finl;
  throw;
}
 
double Champ_Face_PolyMAC_CDO::valeur_a_elem_compo(const DoubleVect& position, int poly, int ncomp) const
{
  Cerr << "Champ_Face_PolyMAC_CDO::" <<__func__ << " is not coded !" << finl;
  throw;
}
 
//interpolation of the integral of v around a dual edge, multiplied by the edge length
void Champ_Face_PolyMAC_CDO::init_ra() const
{
  const Domaine_PolyMAC_CDO& domaine = ref_cast(Domaine_PolyMAC_CDO,domaine_vf());
  const IntTab& a_f = domaine.arete_faces(), &e_f = domaine.elem_faces(), &f_e = domaine.face_voisins();
  const DoubleTab& nf = domaine.face_normales(), &ta = domaine.ta(), &xs = domaine.domaine().coord_sommets(), &xa = dimension < 3 ? xs : domaine.xa(), &xv = domaine.xv();
  const DoubleVect& fs = domaine.face_surfaces(), &ve = domaine.volumes();
  int i, j, k, l, r, e, f, fb, a, idx;
 
  if (radeb.dimension(0)) return;
  init_fcl(), domaine.init_m2(), init_va();
  radeb.resize(1, 2), racf.resize(0, 3);
 
 
  std::map<int, double> rami;
  std::map<int, std::array<double, 3> > ramf;
  for (a = 0; a < xa.dimension_tot(0); radeb.append_line(raji.dimension(0), rajf.dimension(0)), rami.clear(), ramf.clear(), a++)
    {
      //contribution from each face surrounding the edge
      for (i = 0; i < a_f.dimension(1) && (f = a_f(a, i)) >= 0; i++)
        {
          //orientation relative to the edge
          std::array<double, 3> taz = {{ 0, 0, 1 }}, vec; //tangent vector to the edge and its cross product with xv - xa
          vec = domaine.cross(3, dimension, dimension < 3 ? &taz[0] : &ta(a, 0), &xv(f, 0), nullptr, dimension < 3 ? &xs(a, 0): &xa(a, 0));
          int sgn = domaine.dot(&vec[0], &nf(f, 0)) > 0 ? 1 : -1;
 
          //m2 contribution
          for (j = 0; j < 2 && (e = f_e(f, j)) >= 0; j++)
            for (k = domaine.m2d(e), idx = 0; k < domaine.m2d(e + 1); k++, idx++)
              for (l = domaine.m2i(k); f == e_f(e, idx) && l < domaine.m2i(k + 1); l++)
                if (fcl_(fb = e_f(e, domaine.m2j(l)), 0) < 2) rami[fb] += sgn * (e == f_e(f, 0) ? 1 : -1) * (e == f_e(fb, 0) ? 1 : -1) * ve(e) * domaine.m2c(l) / fs(f);
                else if (fcl_(fb, 0) == 3)
                  for (r = 0; r < dimension; r++)
                    ramf[fb][r] += sgn * (e == f_e(f, 0) ? 1 : -1) * (e == f_e(fb, 0) ? 1 : -1) * ve(e) * domaine.m2c(l) / fs(f) * nf(fb, r) / fs(fb);
          //boundary part -> with va if Neumann or Symetrie, with val_imp if Dirichlet_homogene
          if (fcl_(f, 0) == 1 || fcl_(f, 0) == 2)
            {
              for (k = vadeb(a, 0); k < vadeb(a + 1, 0); k++) rami[vaji(k)] += sgn * domaine.dot(&xa(a, 0), &vaci(k, 0), &xv(f, 0));
              for (k = vadeb(a, 1); k < vadeb(a + 1, 1); k++) ramf[vajf(k, 0)][vajf(k, 1)] += sgn * domaine.dot(&xa(a, 0), &vacf(k, 0), &xv(f, 0));
            }
          else if (fcl_(f, 0) == 3)
            for (k = 0; k < dimension; k++)
              ramf[f][k] += sgn * (xa(a, k) - xv(f, k));
        }
      //filling
      double la = dimension < 3 ? 1 : domaine.longueur_aretes()(a);
      for (auto &&kv : rami)
        if (std::fabs(kv.second) > 1e-8 * la)
          raji.append_line(kv.first), raci.append_line(kv.second * la);
      for (auto &&kv : ramf)
        if (domaine.dot(&kv.second[0], &kv.second[0]) > 1e-16 * la * la)
          rajf.append_line(kv.first), racf.append_line(kv.second[0] * la, kv.second[1] * la, kv.second[2] * la);
    }
  CRIMP(radeb), CRIMP(raji), CRIMP(rajf), CRIMP(raci), CRIMP(racf);
}
 
//interpolation of velocity at edges (in the plane normal to each edge)
void Champ_Face_PolyMAC_CDO::init_va() const
{
  const Domaine_PolyMAC_CDO& domaine = ref_cast(Domaine_PolyMAC_CDO,domaine_vf());
  const IntTab& f_e = domaine.face_voisins(), &a_f = domaine.arete_faces();
  const DoubleTab& xv = domaine.xv(), &ta = domaine.ta(), &xs = domaine.domaine().coord_sommets(),
                   &xa = dimension < 3 ? xs : domaine.xa(), &nf = domaine.face_normales(), &xp = domaine.xp();
  const DoubleVect& fs = domaine.face_surfaces(), &pf = equation().milieu().porosite_face(), &pe = equation().milieu().porosite_elem();
 
  int i, j, k, l, m, e, f, fb, a;
 
  if (vadeb.dimension(0)) return;
  domaine.init_m1(), domaine.init_ve();
  std::map<int, std::array<double, 3>> vami;
  std::map<std::array<int, 2>, std::array<double, 3>> vamf;
  vadeb.resize(1, 2), vajf.resize(0, 2), vaci.resize(0, 3), vacf.resize(0, 3);
 
 
  Matrice33 M, iM;
  for (a = 0; a < xa.dimension_tot(0); vadeb.append_line(vaji.dimension(0), vajf.dimension(0)), vami.clear(), vamf.clear(), a++)
    {
      /* list of faces ordered in trigonometric (counter-clockwise) order */
      //local basis
      std::array<double, 3> taz = {{ 0, 0, 1 }}, v0 = domaine.cross(3, dimension, dimension < 3 ? &taz[0] : &ta(a, 0), &xv(a_f(a, 0), 0), nullptr, &xa(a, 0)),
      v1 = domaine.cross(3, 3, dimension < 3 ? &taz[0] : &ta(a, 0), &v0[0]);
      // list of faces, ordered by rotating in the local basis
      std::map<double, int> fmap;
      for (i = 0; i < a_f.dimension(1) && (f = a_f(a, i)) >= 0; i++) fmap[atan2(domaine.dot(&xv(f, 0), &v1[0], &xa(a, 0)), domaine.dot(&xv(f, 0), &v0[0], &xa(a, 0)))] = f;
      std::vector<int> fas, sgn;
      for (auto && kv : fmap) fas.push_back(kv.second);
      //face orientations
      /* compute the interpolation */
      DoubleTrav xe((int)fas.size(), dimension), surf((int)fas.size(), 2), nsurf((int)fas.size(), 2, 3); //intermediate points between xv, surfaces and normals of the dual edge triangles
      double surf_tot = 0, xg[3] = { 0, 0, 0 }; //total surface, center of gravity, normal to the edge
      //orientations of faces relative to the contour
      for (i = 0; i < (int) fas.size(); i++)
        {
          std::array <double, 3> vec = domaine.cross(3, dimension, dimension < 3 ? &taz[0] : &ta(a, 0), &xv(fas[i], 0), nullptr, &xa(a, 0));
          sgn.push_back(domaine.dot(&vec[0], &nf(fas[i], 0)) > 0 ? 1 : -1);
        }
      for (i = 0; i < (int) fas.size(); i++)
        {
          //intermediate points between faces (not necessarily xp)
          int fa[2] = { f = fas[i], fb = fas[(i + 1) * (i + 1 < (int) fas.size())] }, sgfa[2] = { sgn[i], sgn[(i + 1) * (i + 1 < (int) fas.size())] };
          if (domaine.dot(&xv(fb, 0), &nf(f, 0), &xv(f, 0)) * sgn[i] < 0) //angle > 180 -> go through xa
            for (k = 0; k < dimension; k++) xe(i, k) = xa(a, k);
          else if (domaine.dot(&xv(f, 0), &xv(fb, 0), &xa(a, 0), &xa(a, 0)) - (dimension < 3 ? 0 : domaine.dot(&xv(f, 0), &ta(a, 0), &xa(a, 0)) * domaine.dot(&xv(fb, 0), &ta(a, 0), &xa(a, 0))) <= 0) //angle obtus -> parallelogramme
            for (k = 0; k < dimension; k++) xe(i, k) = xv(f, k) + xv(fb, k) - xa(a, k);
          else //acute angle -> intersection of face normals. Solve (xe - xv).(xv - xa) = 0 for both faces, and xe.ta = ta.(xf + xfb) / 2 in 3D
            {
              double s[3];
              for (j = 0; j < 2; j++)
                for (k = 0, s[j] = domaine.dot(&xv(fa[j], 0), &xv(fa[j], 0), &xa(a, 0), &xa(a, 0)); k < 3; k++)
                  M(j, k) = k < dimension ? xv(fa[j], k) - xa(a, k) : 0;
              for (k = 0, s[2] = dimension < 3 ? 0 : (domaine.dot(&xv(f, 0), &ta(a, 0), &xa(a, 0)) + domaine.dot(&xv(fb, 0), &ta(a, 0), &xa(a, 0))) / 2; k < 3; k++)
                M(2, k) = dimension < 3 ? (k == 2) : ta(a, k);
              double eps = Matrice33::inverse(M, iM, 0);
              for (j = 0; eps != 0 && j < dimension; j++)
                for (k = 0, xe(i, j) = xa(a, j); k < dimension; k++) xe(i, j) += iM(j, k) * s[k];
              //if xe ends up on the wrong side of a face, use xp (if the element exists) or xf + xfb / 2 (otherwise)
              if (eps == 0 || domaine.dot(&xe(i, 0), &nf(f, 0), &xv(f, 0)) * sgfa[0] < 0 || domaine.dot(&xe(i, 0), &nf(fb, 0), &xv(fb, 0)) * sgfa[1] > 0)
                {
                  if ((e = f_e(f, sgfa[0] > 0)) >= 0)
                    for (k = 0; k < dimension; k++) xe(i, k) = xp(e, k);
                  else for (k = 0; k < dimension; k++) xe(i, k) = (xv(f, k) + xv(fb, k)) / 2;
                }
            }
 
          //surfaces and center of gravity
          for (j = 0; j < 2; j++)
            {
              std::array<double, 3> vsurf3 = domaine.cross(dimension, dimension, &xe(i, 0), &xv(fa[j], 0), &xa(a, 0), &xa(a, 0));
              //orientation relative to ta
              for (k = 0, l = (dimension < 3 ? vsurf3[2] < 0 : domaine.dot(&ta(a, 0), &vsurf3[0]) < 0); k < 3; k++) vsurf3[k] *=  (l ? -1. : 1.) / 2;
              surf(i, j) = dimension < 3 ? std::fabs(vsurf3[2]) : sqrt(domaine.dot(&vsurf3[0], &vsurf3[0]));
              surf_tot += surf(i, j);
              for (k = 0; surf(i, j) > 1e-16 && k < 3; k++) nsurf(i, j, k) = vsurf3[k] / surf(i, j);
              for (k = 0; k < dimension; k++) xg[k] += surf(i, j) * (xv(fa[j], k) + xe(i, k) - 2 * xa(a, k)) / 3;
            }
        }
      for (k = 0; k < dimension; k++) xg[k] = xa(a, k) + (surf_tot ==0 ? 0. : xg[k] / surf_tot);
 
      if (!surf_tot) continue;
      M = Matrice33(surf_tot, 0, 0, 0, surf_tot, 0, 0, 0, surf_tot);
      /* loops over the triangles of the dual edge */
      std::array<double, 3> vec;
      for (i = 0; i < (int) fas.size(); i++)
        for (j = 0; j < 2; j++)
          if (surf(i, j) > 1e-16)
            {
              f = fas[(i + j) * (i + j < (int) fas.size())], e = f_e(f, sgn[i] > 0); //face, element (if it exists)
              double x[2][3] = { { 0, }, }; //boundary points: xv(fa[i]) -> xe(i) (j = 0) or xe(i) -> xv(fa[i + 1]) (j = 1)
              for (k = 0; k < 2; k++)
                for (l = 0; l < dimension; l++) x[k][l] = j == k ? xv(f, l) : xe(i, l);
 
              //edge-point segments: using va
              for (k = 0; k < 2; k++)
                {
                  for (l = 0; l < dimension; l++) vec[l] = (xa(a, l) + x[k][l]) / 2 - xg[l];
                  vec = domaine.cross(3, dimension, &nsurf(i, j, 0), &vec[0]);
                  for (l = 0; l < 3; l++)
                    for (m = 0; m < dimension; m++) M(l, m) -= (k ? -1 : 1) * (x[k][m] - xa(a, m)) * vec[l];
                }
 
              //segment x[0] -> x[1]
              for (k = 0; k < dimension; k++) vec[k] = (x[0][k] + x[1][k]) / 2 - xg[k];
              vec = domaine.cross(3, dimension, &nsurf(i, j, 0), &vec[0]);
              //vf part
              if (fcl_(f, 0) < 2)
                for (k = 0; k < dimension; k++)
                  vami[f][k] += domaine.dot(x[1], &nf(f, 0), x[0]) / fs(f) * vec[k];
              else if (fcl_(f, 0) == 3)
                for (k = 0; k < dimension; k++)
                  for (l = 0; l < dimension; l++)
                    vamf[ {{ f, l }}][l] += domaine.dot(x[1], &nf(f, 0), x[0]) / fs(f) * vec[k] * nf(f, l) / fs(f);
              //complement: using ve (e >= 0) / imposed value (Dirichlet) / va (Neumann/Symmetry)
              if (e >= 0)
                for (k = domaine.vedeb(e); k < domaine.vedeb(e + 1); k++) //ve
                  {
                    if (fcl_(fb = domaine.veji(k), 0) < 2)
                      for (l = 0; l < dimension; l++)
                        vami[fb][l] += (domaine.dot(x[1], &domaine.veci(k, 0), x[0]) - domaine.dot(x[1], &nf(f, 0), x[0]) * domaine.dot(&domaine.veci(k, 0), &nf(f, 0)) / (fs(f) * fs(f))) * vec[l] * pf(fb) / pe(e);
                    else if (fcl_(fb, 0) == 3)
                      for (l = 0; l < dimension; l++)
                        for (m = 0; m < dimension; m++)
                          vamf[ {{ fb, m }}][l] += (domaine.dot(x[1], &domaine.veci(k, 0), x[0]) - domaine.dot(x[1], &nf(f, 0), x[0]) * domaine.dot(&domaine.veci(k, 0), &nf(f, 0)) / (fs(f) * fs(f)))
                    * vec[l] * nf(fb, m) / fs(f) * pf(fb) / pe(e);
                  }
              else if (fcl_(f, 0) == 3)
                for (k = 0; k < dimension; k++)
                  for (l = 0; l < dimension; l++) //imposed value
                    vamf[ {{ f, k }}][l] += (x[1][k] - x[0][k] - domaine.dot(x[1], &nf(f, 0), x[0]) * nf(f, k) / (fs(f) * fs(f))) * vec[l];
              else if (fcl_(f, 0) == 1 || fcl_(f, 0) == 2)
                for (k = 0; k < dimension; k++)
                  for (l = 0; l < dimension; l++) //va
                    M(k, l) -= (x[1][l] - x[0][l] - domaine.dot(x[1], &nf(f, 0), x[0]) * nf(f, l) / (fs(f) * fs(f))) * vec[k];
 
              //normal part of the triangle (3D only): also using vf/ve/imposed value
              if (dimension < 3) continue;
              else if (e >= 0)
                for (k = domaine.vedeb(e); k < domaine.vedeb(e + 1); k++) //ve
                  {
                    if (fcl_(fb = domaine.veji(k), 0) < 2)
                      for (l = 0; l < dimension; l++)
                        vami[fb][l] += surf(i, j) * domaine.dot(&nsurf(i, j, 0), &domaine.veci(k, 0)) * nsurf(i, j, l) * pf(fb) / pe(e);
                    else if (fcl_(fb, 0) == 3)
                      for (l = 0; l < dimension; l++)
                        for (m = 0; m < dimension; m++)
                          vamf[ {{ fb, m }}][l] += surf(i, j) * domaine.dot(&nsurf(i, j, 0), &domaine.veci(k, 0)) * nsurf(i, j, l) * nf(fb, m) / fs(f) * pf(fb) / pe(e);
                  }
              else if (fcl_(f, 0) == 3)
                for (k = 0; k < dimension; k++)
                  for (l = 0; l < dimension; l++) //imposed value
                    vamf[ {{ f, k }}][l] += surf(i, j) * nsurf(i, j, k) * nsurf(i, j, l);
              else if (fcl_(f, 0) == 1 || fcl_(f, 0) == 2)
                for (k = 0; k < dimension; k++)
                  for (l = 0; l < dimension; l++) //va
                    M(k, l) -= surf(i, j) * nsurf(i, j, k) * nsurf(i, j, l);
            }
 
      //inversion of M and filling
      double eps = Matrice33::inverse(M, iM, 0);
      if (std::fabs(eps) < 1e-24) continue;
      for (auto && kv : vami)
        {
          double inc[3] = { 0, 0, 0 };
          for (k = 0; k < dimension; k++)
            for (l = 0; l < dimension; l++) inc[k] += iM(k, l) * kv.second[l];
          if (domaine.dot(inc, inc) > 1e-16) vaji.append_line(kv.first), vaci.append_line(inc[0], inc[1], inc[2]);
        }
      for (auto && kv : vamf)
        {
          double inc[3] = { 0, 0, 0 };
          for (k = 0; k < dimension; k++)
            for (l = 0; l < dimension; l++) inc[k] += iM(k, l) * kv.second[l];
          if (domaine.dot(inc, inc) > 1e-16) vajf.append_line(kv.first[0], kv.first[1]), vacf.append_line(inc[0], inc[1], inc[2]);
        }
    }
  CRIMP(vadeb), CRIMP(vaji), CRIMP(vajf), CRIMP(vaci), CRIMP(vacf);
}
 
/* velocity at elements */
void Champ_Face_PolyMAC_CDO::interp_ve(const DoubleTab& inco, DoubleTab& val, bool is_vit) const
{
  const Domaine_PolyMAC_CDO& domaine = ref_cast(Domaine_PolyMAC_CDO,domaine_vf());
  int e, f, j, k, r;
 
  domaine.init_ve();
  val = 0;
 
  if (polymac_flica5)
    {
      const DoubleVect& pf = equation().milieu().porosite_face(), &pe = equation().milieu().porosite_elem();
      for (e = 0; e < val.dimension(0); e++)
        for (j = domaine.vedeb(e); j < domaine.vedeb(e + 1); j++)
          {
            f = domaine.veji(j);
            const double coef = is_vit ? pf(f) / pe(e) : 1.0;
            for (r = 0; r < dimension; r++) val(e, r) += domaine.veci(j, r) * inco(f) * coef;
          }
    }
  else
    {
      const Conds_lim& cls = domaine_Cl_dis().les_conditions_limites();
      const DoubleTab& nf = domaine.face_normales();
      const DoubleVect& fs = domaine.face_surfaces(), *pf = mon_equation_non_nul() ? &equation().milieu().porosite_face() : nullptr, *pe = pf ? &equation().milieu().porosite_elem() : nullptr;
 
      for (e = 0; e < val.dimension(0); e++)
        for (j = domaine.vedeb(e); j < domaine.vedeb(e + 1); j++)
          if (fcl_(f = domaine.veji(j), 0) < 2) //computed velocity
            {
              const double coef = is_vit && pf ? (*pf)(f) / (*pe)(e) : 1.0;
              for (r = 0; r < dimension; r++) val(e, r) += domaine.veci(j, r) * inco(f) * coef;
            }
          else if (fcl_(f, 0) == 3)
            for (k = 0; k < dimension; k++)
              for (r = 0; r < dimension; r++) //Dirichlet
                {
                  const double coef = is_vit && pf ? (*pf)(f) / (*pe)(e) : 1.0;
                  val(e, r) += domaine.veci(j, r) * ref_cast(Dirichlet, cls[fcl_(f, 1)].valeur()).val_imp(fcl_(f, 2), k) * nf(f, k) / fs(f) * coef;
                }
    }
}
 
/* velocity at elements over a list of elements */
void Champ_Face_PolyMAC_CDO::interp_ve(const DoubleTab& inco, const IntVect& les_polys, DoubleTab& val, bool is_vit) const
{
  const Domaine_PolyMAC_CDO& domaine = ref_cast(Domaine_PolyMAC_CDO,domaine_vf());
  int e, f, j, r;
 
  domaine.init_ve();
 
  const DoubleVect *pf = mon_equation_non_nul() ? &equation().milieu().porosite_face() : nullptr, *pe = pf ? &equation().milieu().porosite_elem() : nullptr;
  for (int poly = 0; poly < les_polys.size(); poly++)
    {
      e = les_polys(poly);
      if (e!=-1)
        {
          for (r = 0; r < dimension; r++) val(e, r) = 0;
          for (j = domaine.vedeb(e); j < domaine.vedeb(e + 1); j++)
            {
              f = domaine.veji(j);
              const double coef = is_vit && pf ? (*pf)(f) / (*pe)(e) : 1.0;
              for (r = 0; r < dimension; r++) val(e, r) += domaine.veci(j, r) * inco(f) * coef;
            }
        }
    }
}
 
/* gradient d_j v_i at elements */
void Champ_Face_PolyMAC_CDO::interp_gve(const DoubleTab& inco, DoubleTab& vals) const
{
  const Domaine_PolyMAC_CDO& domaine = ref_cast(Domaine_PolyMAC_CDO,domaine_vf());
  const Conds_lim& cls = domaine_Cl_dis().les_conditions_limites();
  const IntTab& f_e = domaine.face_voisins();
  const DoubleTab& nf = domaine.face_normales();
  const DoubleVect& fs = domaine.face_surfaces();
  int i, j, k, f, e;
  double scal;
 
  domaine.init_m2solv();
  //velocities at elements, gradient of each component at faces
  DoubleTrav ve1(0, dimension);
  domaine.domaine().creer_tableau_elements(ve1);
  interp_ve(inco, ve1);
 
  //gradient at dual faces, field (nf.grad)v (obtained by inverting m2)
  DoubleTrav gv(domaine.nb_faces_tot(), dimension), cgv, ngv;
  domaine.creer_tableau_faces(cgv), domaine.creer_tableau_faces(ngv);
  for (f = 0; f < domaine.nb_faces_tot(); f++)
    if (fcl_(f, 0) != 1) //gve = 0 if Neumann
      {
        for (i = 0; i < 2; i++)
          if ((e = f_e(f, i)) >= 0)
            for (k = 0; k < dimension; k++) gv(f, k) += (i ? 1 : -1) * fs(f) * ve1(e, k); //internal element -> using ve1
          else if (fcl_(f, 0) == 3)
            for (k = 0; k < dimension; k++) //Dirichlet boundary -> using imposed value
              gv(f, k) += (i ? 1 : -1) * fs(f) * ref_cast(Dirichlet, cls[fcl_(f, 1)].valeur()).val_imp(fcl_(f, 2), k);
        //if Symmetry -> keep only the normal component at the face
        if (fcl_(f, 0) == 2)
          for (k = 0, scal = domaine.dot(&nf(f, 0), &gv(f, 0)) / fs(f); k < dimension; k++) gv(f, k) = scal * nf(f, k) / fs(f);
      }
 
  //for each component:
  for (k = 0; k < dimension; k++)
    {
      /* invert m2 -> to obtain (nf.grad) v_k */
      for (f = 0; f < domaine.nb_faces_tot(); f++) cgv(f) = gv(f, k);
      domaine.m2solv.resoudre_systeme(domaine.m2mat, cgv, ngv);
      /* re-interpolation at elements -> grad v_k */
      for (e = 0; e < domaine.nb_elem(); e++)
        for (i = domaine.vedeb(e); i < domaine.vedeb(e + 1); i++)
          for (j = 0; j < dimension; j++)
            vals.addr()[(e * dimension + k) * dimension + j] += domaine.veci(i, j) * ngv(domaine.veji(i));
    }
 
  vals.echange_espace_virtuel();
}
 
DoubleTab& Champ_Face_PolyMAC_CDO::valeur_aux_elems_(const DoubleTab& val_face, const DoubleTab& positions, const IntVect& les_polys, DoubleTab& val_elem) const
{
  const Champ_base& cha = le_champ();
  int nb_compo = cha.nb_comp(), N = val_face.line_size(), D = dimension;
  assert(val_elem.line_size() == nb_compo * N);
  assert((positions.dimension(0) == les_polys.size()) || (positions.dimension_tot(0) == les_polys.size()));
 
  if (val_elem.nb_dim() > 2)
    Process::exit("TRUST error in Champ_Face_PolyMAC_CDO::valeur_aux_elems_ : The DoubleTab val has more than 2 entries !");
 
  if (nb_compo == 1)
    Process::exit("TRUST error in Champ_Face_PolyMAC_CDO::valeur_aux_elems_ : A scalar field cannot be of Champ_Face type !");
 
  // only if Champ_Face_PolyMAC_CDO because interp_ve needs mon_dom_cl_dis ...
  if (!mon_dom_cl_dis && que_suis_je() == "Champ_Face_PolyMAC_CDO")
    return val_elem; //nothing can be done until the boundary conditions are known
 
  //interpolate ve over all elements, then restrict to les_polys
  DoubleTrav ve(0, N * D);
  const Domaine_VF& domdom = ref_cast(Domaine_VF, domaine_vf());
  domdom.domaine().creer_tableau_elements(ve);
 
  if (polymac_flica5 && que_suis_je() == "Champ_Face_PolyMAC_CDO")
    {
      bool is_vit = false && cha.le_nom().debute_par("vitesse");
      interp_ve(val_face, les_polys, ve, is_vit);
    }
  else
    {
      bool is_vit = cha.le_nom().debute_par("vitesse") && !cha.le_nom().debute_par("vitesse_debitante");
      interp_ve(val_face, ve, is_vit);
    }
 
  for (int p = 0; p < les_polys.size(); p++)
    for (int r = 0, e = les_polys(p); e < domdom.nb_elem() && r < N * D; r++)
      val_elem(p, r) = (e == -1) ? 0. : ve(e, r);
  return val_elem;
}
 
DoubleTab& Champ_Face_PolyMAC_CDO::valeur_aux_elems(const DoubleTab& positions, const IntVect& les_polys, DoubleTab& val_elem) const
{
  return valeur_aux_elems_(le_champ().valeurs(), positions, les_polys, val_elem);
}
 
DoubleTab& Champ_Face_PolyMAC_CDO::valeur_aux_elems_passe(const DoubleTab& positions, const IntVect& les_polys, DoubleTab& val_elem) const
{
  return valeur_aux_elems_(le_champ().passe(), positions, les_polys, val_elem);
}
 
DoubleVect& Champ_Face_PolyMAC_CDO::valeur_aux_elems_compo(const DoubleTab& positions, const IntVect& les_polys, DoubleVect& val, int ncomp) const
{
  if (!polymac_flica5 || que_suis_je() != "Champ_Face_PolyMAC_CDO") fcl();
  const Champ_base& cha=le_champ();
  assert(val.size_totale() >= les_polys.size());
 
  // only if Champ_Face_PolyMAC_CDO because interp_ve needs mon_dom_cl_dis ...
  if (!mon_dom_cl_dis && que_suis_je() ==  "Champ_Face_PolyMAC_CDO")
    return val; //nothing can be done until the boundary conditions are known
 
  //interpolate ve over all elements, then restrict to les_polys
  DoubleTrav ve(0, dimension * cha.valeurs().line_size());
  ref_cast(Domaine_VF,domaine_vf()).domaine().creer_tableau_elements(ve);
  if (polymac_flica5 && que_suis_je() == "Champ_Face_PolyMAC_CDO")
    {
      bool is_vit = false && cha.le_nom().debute_par("vitesse");
      interp_ve(cha.valeurs(), les_polys, ve, is_vit);
    }
  else
    interp_ve(cha.valeurs(), ve);
 
  for (int p = 0; p < les_polys.size(); p++) val(p) = (les_polys(p) == -1) ? 0. : ve(les_polys(p), ncomp);
 
  return val;
}
 
DoubleTab& Champ_Face_PolyMAC_CDO::remplir_coord_noeuds(DoubleTab& positions) const
{
  Cerr << "Champ_Face_PolyMAC_CDO::" <<__func__ << " is not coded !" << finl;
  throw;
}
 
int Champ_Face_PolyMAC_CDO::remplir_coord_noeuds_et_polys(DoubleTab& positions, IntVect& polys) const
{
  Cerr << "Champ_Face_PolyMAC_CDO::" <<__func__ << " is not coded !" << finl;
  throw;
}
 
DoubleTab& Champ_Face_PolyMAC_CDO::valeur_aux_faces(DoubleTab& val) const
{
  const Champ_base& cha=le_champ();
  int nb_compo=cha.nb_comp();
 
  if (nb_compo == 1)
    Process::exit("Champ_Face_PolyMAC_CDO::valeur_aux_faces : A scalar field cannot be of Champ_Face type !");
 
  const Domaine_PolyMAC_CDO& domaine = ref_cast(Domaine_PolyMAC_CDO,domaine_vf());
  val.resize(domaine.nb_faces(), dimension), val = 0;
 
  for (int f = 0; f < domaine.nb_faces(); f++)
    for (int r = 0; r < dimension; r++) val(f, r) = cha.valeurs()(f) * domaine.face_normales(f, r) / domaine.face_surfaces(f);
  return val;
}
 
DoubleVect& Champ_Face_PolyMAC_CDO::calcul_S_barre_sans_contrib_paroi(const DoubleTab& vitesse, DoubleVect& SMA_barre) const
{
  // with boundary contribution for now...
  abort();
  return calcul_S_barre(vitesse, SMA_barre);
}
 
DoubleVect& Champ_Face_PolyMAC_CDO::calcul_S_barre(const DoubleTab& vitesse, DoubleVect& SMA_barre) const
{
  const Domaine_PolyMAC_CDO& domaine = ref_cast(Domaine_PolyMAC_CDO, domaine_vf());
  DoubleTab duidxj(0, dimension, dimension);
  domaine.domaine().creer_tableau_elements(duidxj);
  interp_gve(vitesse, duidxj);
 
  for (int elem = 0; elem < domaine.domaine().nb_elem(); elem++)
    {
      double temp = 0.;
      for (int i = 0; i < dimension; i++)
        for (int j = 0; j < dimension; j++)
          {
            double Sij = 0.5 * (duidxj(elem, i, j) + duidxj(elem, j, i));
            temp += Sij * Sij;
          }
      SMA_barre(elem) = 2. * temp;
    }
 
  return SMA_barre;
 
}
 
DoubleTab& Champ_Face_PolyMAC_CDO::trace(const Frontiere_dis_base& fr, DoubleTab& x, double t, int distant) const
{
  assert(distant==0);
  init_fcl();
  const bool vectoriel = (le_champ().nb_comp() > 1);
  const int dim = vectoriel ? dimension : 1;
  const Front_VF& fr_vf = ref_cast(Front_VF, fr);
  const Domaine_PolyMAC_CDO& domaine = ref_cast(Domaine_PolyMAC_CDO,domaine_vf());
  const IntTab& face_voisins = domaine.face_voisins();
  const DoubleTab& val = valeurs(t);
  DoubleTrav ve(0, dimension);
  if (vectoriel)
    {
      domaine.domaine().creer_tableau_elements(ve);
      interp_ve(val, ve);
    }
 
  for (int i = 0; i < fr_vf.nb_faces(); i++)
    {
      const int face = fr_vf.num_premiere_face() + i;
      for (int dir = 0; dir < 2; dir++)
        {
          const int elem = face_voisins(face, dir);
          if (elem != -1)
            {
              for (int d = 0; d < dim; d++) x(i, d) = vectoriel ? ve(elem, d) : val[face];
            }
        }
    }
 
  return x;
}
 
int Champ_Face_PolyMAC_CDO::nb_valeurs_nodales() const
{
  ConstDoubleTab_parts val_part(valeurs());
  return val_part[0].dimension(0) + val_part[1].dimension(0);
}