/*============================================================================
*
*                    Code_Saturne version 1.3
*                    ------------------------
*
*
*     This file is part of the Code_Saturne Kernel, element of the
*     Code_Saturne CFD tool.
*
*     Copyright (C) 1998-2007 EDF S.A., France
*
*     contact: saturne-support@edf.fr
*
*     The Code_Saturne Kernel is free software; you can redistribute it
*     and/or modify it under the terms of the GNU General Public License
*     as published by the Free Software Foundation; either version 2 of
*     the License, or (at your option) any later version.
*
*     The Code_Saturne Kernel is distributed in the hope that it will be
*     useful, but WITHOUT ANY WARRANTY; without even the implied warranty
*     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*     GNU General Public License for more details.
*
*     You should have received a copy of the GNU General Public License
*     along with the Code_Saturne Kernel; if not, write to the
*     Free Software Foundation, Inc.,
*     51 Franklin St, Fifth Floor,
*     Boston, MA  02110-1301  USA
*
*============================================================================*/

/*============================================================================
 * Compute several mesh quality criteria.
 *============================================================================*/

/*----------------------------------------------------------------------------
 * Standard C library headers
 *----------------------------------------------------------------------------*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>

/*----------------------------------------------------------------------------
 * BFT library headers
 *----------------------------------------------------------------------------*/

#include <bft_mem.h>
#include <bft_error.h>
#include <bft_printf.h>

/*----------------------------------------------------------------------------
 * Local headers
 *----------------------------------------------------------------------------*/

#include "cs_maillage.h"
#include "cs_maillage_grd.h"
#include "cs_post.h"

/*----------------------------------------------------------------------------
 * Header for the current file
 *----------------------------------------------------------------------------*/

#include "cs_mesh_quality.h"

/*----------------------------------------------------------------------------*/

#ifdef __cplusplus
extern "C" {
#if 0
} /* Fake brace to force Emacs auto-indentation back to column 0 */
#endif
#endif /* __cplusplus */

/*=============================================================================
 * Local Macro Definitions
 *============================================================================*/

#define CS_MESH_QUALITY_N_SUBS  10

#undef _CROSS_PRODUCT_3D
#undef _DOT_PRODUCT_3D
#undef _MODULE_3D
#undef _COSINE_3D

#define _CROSS_PRODUCT_3D(cross_v1_v2, v1, v2) ( \
 cross_v1_v2[0] = v1[1]*v2[2] - v1[2]*v2[1],   \
 cross_v1_v2[1] = v1[2]*v2[0] - v1[0]*v2[2],   \
 cross_v1_v2[2] = v1[0]*v2[1] - v1[1]*v2[0]  )

#define _DOT_PRODUCT_3D(v1, v2) ( \
 v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2])

#define _MODULE_3D(v) \
 sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2])

#define _COSINE_3D(v1, v2) (\
 _DOT_PRODUCT_3D(v1, v2) / (_MODULE_3D(v1) * _MODULE_3D(v2)) )

/*============================================================================
 * Static global variables
 *============================================================================*/

/*============================================================================
 * Private function definitions
 *============================================================================*/

/*----------------------------------------------------------------------------
 * Compute the minimum and the maximum of a vector (locally).
 *
 * parameters:
 *   n_vals    --> local number of elements
 *   var       --> pointer to vector
 *   min       <-- minimum
 *   max       <-- maximum
 *----------------------------------------------------------------------------*/

static void
_compute_local_minmax(cs_int_t         n_vals,
                      const cs_real_t  var[],
                      cs_real_t       *min,
                      cs_real_t       *max)
{
  cs_int_t  i;
  cs_real_t  _min = var[0], _max = var[0];

  for (i = 1; i < n_vals; i++) {
    _min = CS_MIN(_min, var[i]);
    _max = CS_MAX(_max, var[i]);
  }

  if (min != NULL)  *min = _min;
  if (max != NULL)  *max = _max;
}

/*----------------------------------------------------------------------------
 * Display the distribution of values of a real vector.
 *
 * parameters:
 *   n_vals    --> local number of elements
 *   var       --> pointer to vector
 *----------------------------------------------------------------------------*/

static void
_display_histograms(cs_int_t         n_vals,
                    const cs_real_t  var[])
{
  cs_int_t  i_step, i_val, counter;

  cs_real_t  step;
  cs_real_t  l_bound, u_bound;
  cs_real_t  max, min;

  cs_real_t  _max = var[0], _min = var[0];

  const cs_int_t  n_steps = CS_MESH_QUALITY_N_SUBS;

  assert (sizeof (double) == sizeof (cs_real_t));

  /* Compute global min and max */

  _compute_local_minmax(n_vals, var, &_min, &_max);

  /* Default initialization */

  min = _min;
  max = _max;

#if defined(_CS_HAVE_MPI)

  if (cs_glob_base_nbr > 1) {
    MPI_Allreduce(&_min, &min, 1, CS_MPI_REAL, MPI_MIN,
                  cs_glob_base_mpi_comm);

    MPI_Allreduce(&_max, &max, 1, CS_MPI_REAL, MPI_MAX,
                  cs_glob_base_mpi_comm);
  }

#endif

  bft_printf(_("    valeur minimale =       %10.5e\n"), min);
  bft_printf(_("    valeur maximale =       %10.5e\n"), max);

  /* Define axis subdivisions */

  if (CS_ABS(max - min) > 0) {

    step = CS_ABS(max - min) / n_steps;
    bft_printf(_("    valeur de l'incrément = %10.5e\n\n"), step);

    /* Number of elements in each subdivision */

    u_bound = min;

    for (i_step = 1; i_step < n_steps; i_step++) {

      l_bound = u_bound;
      u_bound = l_bound + step;

      counter = 0;

      for (i_val = 0; i_val < n_vals; i_val++)
        if (var[i_val] >= l_bound && var[i_val] < u_bound)  counter++;

      bft_printf("    %3d : [ %10.5e ; %10.5e [ = %10d\n",
                 i_step, l_bound, u_bound, counter);

    }

    l_bound = max - step;
    u_bound = max;

    counter = 0;

    for (i_val = 0; i_val < n_vals; i_val++)
      if (var[i_val] >= l_bound)  counter++;

    bft_printf("    %3d : [ %10.5e ; %10.5e [ = %10d\n",
               n_steps, l_bound, u_bound, counter);
  }

}

/*----------------------------------------------------------------------------
 * Compute weighting coefficient and center offsetting coefficient
 * for internal faces.
 *
 * parameters:
 *   mesh             --> pointer to mesh structure.
 *   mesh_quantities  --> pointer to mesh quantities structures.
 *   weighting        <-> array for weigthing coefficient.
 *   offsetting       <-> array for offsetting coefficient.
 *----------------------------------------------------------------------------*/

static void
_compute_weighting_offsetting(const cs_maillage_t      *mesh,
                              const cs_maillage_grd_t  *mesh_quantities,
                              cs_real_t                 weighting[],
                              cs_real_t                 offsetting[])
{
  cs_int_t  i_fac, i_dim;
  cs_int_t  cell1, cell2;

  cs_real_t  intersection;
  cs_real_t  cell_center1[3], cell_center2[3];
  cs_real_t  face_center[3], face_normal[3];
  cs_real_t  v0[3], v1[3], v2[3];

  double  coef0 = 0, coef1 = 0, coef2 = 0;

  const cs_int_t  dim = mesh->dim;

  /* Compute weighting coefficient */
  /*-------------------------------*/

  /* Loop on internal faces */

  for (i_fac = 0; i_fac < mesh->nbr_fac; i_fac++) {

    /* Get local number of the cells in contact with the face */

    cell1 = mesh->fac_cel[2 * i_fac] - 1;
    cell2 = mesh->fac_cel[2 * i_fac + 1] - 1;

    /* Get information on mesh quantities */

    for (i_dim = 0; i_dim < dim; i_dim++) {

      /* Center of gravity for each cell */

      cell_center1[i_dim] = mesh_quantities->cen_cel[cell1*dim + i_dim];
      cell_center2[i_dim] = mesh_quantities->cen_cel[cell2*dim + i_dim];

      /* Face center coordinates */

      face_center[i_dim] = mesh_quantities->cdg_fac[i_fac*dim + i_dim];

      /* Surface vector (orthogonal to the face) */

      face_normal[i_dim] = mesh_quantities->surf_fac[i_fac*dim + i_dim];

    }

    /* Compute weighting coefficient with two approaches. Keep the max value. */

    for (i_dim = 0; i_dim < dim; i_dim++) {

      v0[i_dim] = cell_center2[i_dim] - cell_center1[i_dim];
      v1[i_dim] = face_center[i_dim] - cell_center1[i_dim];
      v2[i_dim] = cell_center2[i_dim] - face_center[i_dim];

    }

    coef0 = _DOT_PRODUCT_3D(v0, face_normal);
    coef1 = _DOT_PRODUCT_3D(v1, face_normal)/coef0;
    coef2 = _DOT_PRODUCT_3D(v2, face_normal)/coef0;

    weighting[i_fac] = CS_MAX(coef1, coef2);

    /* Compute center offsetting coefficient */
    /*---------------------------------------*/

    /* Compute intersection between face and segment defined by the two cell
       centers */

    for (i_dim = 0; i_dim < dim; i_dim++) {

      intersection =  (1 - weighting[i_fac]) * cell_center1[i_dim]
                    + weighting[i_fac] * cell_center2[i_dim];

      v1[i_dim] = intersection - face_center[i_dim];
      v2[i_dim] = cell_center2[i_dim] - cell_center1[i_dim];

    }

    offsetting[i_fac] = _MODULE_3D(v1) / _MODULE_3D(v2);

  } /* End of loop on faces */

}

/*----------------------------------------------------------------------------
 * Compute angle between face normal and segment based on centers of the
 * adjacent cells. Evaluates a level of non-orthogonality.
 *
 * parameters:
 *   mesh             --> pointer to mesh structure.
 *   mesh_quantities  --> pointer to mesh quantities structures.
 *   fac_ortho        <-> array for internal faces.
 *   fbr_ortho        <-> array for border faces.
 *----------------------------------------------------------------------------*/

static void
_compute_orthogonality(const cs_maillage_t      *mesh,
                       const cs_maillage_grd_t  *mesh_quantities,
                       cs_real_t                 fac_ortho[],
                       cs_real_t                 fbr_ortho[])
{
  cs_int_t  i_fac, i_dim;
  cs_int_t  cell1, cell2;

  double  cos_alpha;
  cs_real_t  cell_center1[3], cell_center2[3];
  cs_real_t  face_center[3];
  cs_real_t  face_normal[3], vect[3];

  const double  rad_to_deg = 180. / acos(-1.);
  const cs_int_t  dim = mesh->dim;

  /* Loop on internal faces */
  /*------------------------*/

  for (i_fac = 0; i_fac < mesh->nbr_fac; i_fac++) {

    /* Get local number of the cells beside the face */

    cell1 = mesh->fac_cel[2 * i_fac] - 1;
    cell2 = mesh->fac_cel[2 * i_fac + 1] - 1;

    /* Get information on mesh quantities */

    for (i_dim = 0; i_dim < dim; i_dim++) {

      /* Center of gravity for each cell */

      cell_center1[i_dim] = mesh_quantities->cen_cel[cell1*dim + i_dim];
      cell_center2[i_dim] = mesh_quantities->cen_cel[cell2*dim + i_dim];

      /* Surface vector (orthogonal to the face) */

      face_normal[i_dim] = mesh_quantities->surf_fac[i_fac*dim + i_dim];

    }

    /* Compute angle which evaluates the non-orthogonality. */

    for (i_dim = 0; i_dim < dim; i_dim++)
      vect[i_dim] = cell_center2[i_dim] - cell_center1[i_dim];

    cos_alpha = _COSINE_3D(vect, face_normal);
    cos_alpha = CS_ABS(cos_alpha);

    fac_ortho[i_fac] = acos(cos_alpha) * rad_to_deg;

  } /* End of loop on internal faces */

  /* Loop on border faces */
  /*----------------------*/

  for (i_fac = 0; i_fac < mesh->nbr_fbr; i_fac++) {

    /* Get local number of the cell beside the face */

    cell1 = mesh->fbr_cel[i_fac] - 1;

    /* Get information on mesh quantities */

    for (i_dim = 0; i_dim < dim; i_dim++) {

      /* Center of gravity of the cell */

      cell_center1[i_dim] = mesh_quantities->cen_cel[cell1*dim + i_dim];

      /* Face center coordinates */

      face_center[i_dim] = mesh_quantities->cdg_fbr[i_fac*dim + i_dim];

      /* Surface vector (orthogonal to the face) */

      face_normal[i_dim] = mesh_quantities->surf_fbr[i_fac*dim + i_dim];

    }

    /* Compute alpha: angle wich evaluate the difference with orthogonality. */

    for (i_dim = 0; i_dim < dim; i_dim++)
      vect[i_dim] = face_center[i_dim] - cell_center1[i_dim];

    cos_alpha = _COSINE_3D(vect, face_normal);
    cos_alpha = CS_ABS(cos_alpha);

    fbr_ortho[i_fac] = acos(cos_alpha) * rad_to_deg;

  } /* End of loop on border faces */

}

/*----------------------------------------------------------------------------
 * Evaluate face warping angle.
 *
 * parameters:
 *   idx_start       --> first vertex index
 *   idx_end         --> last vertex index
 *   face_vertex_num --> face -> vertices connectivity
 *   face_normal     --> face normal
 *   vertex_coords   --> vertices coordinates
 *   face_warping    <-- face warping angle
 *----------------------------------------------------------------------------*/

static void
_get_face_warping(cs_int_t         idx_start,
                  cs_int_t         idx_end,
                  const cs_real_t  face_normal[],
                  const cs_int_t   face_vertex_num[],
                  const cs_real_t  vertex_coords[],
                  double           face_warping[])
{
  cs_int_t  i_dim, idx;
  cs_int_t  vertex_num1, vertex_num2;

  double  edge_cos_alpha;
  cs_real_t  vect[3];

  double  cos_alpha = 1.;

  const int dim = 3;
  const double  rad_to_deg = 180. / acos(-1.);

  /* Loop on edges */

  for (idx = idx_start; idx < idx_end - 1; idx++) {

    vertex_num1 = face_vertex_num[idx] - 1;
    vertex_num2 = face_vertex_num[idx + 1] - 1;

    /* Get vertex coordinates */

    for (i_dim = 0; i_dim < dim; i_dim++)
      vect[i_dim] =  vertex_coords[vertex_num2*dim + i_dim]
                   - vertex_coords[vertex_num1*dim + i_dim];

    edge_cos_alpha = _COSINE_3D(vect, face_normal);
    edge_cos_alpha = CS_ABS(edge_cos_alpha);
    cos_alpha = CS_MIN(cos_alpha, edge_cos_alpha);

  }

  /* Last edge */

  vertex_num1 = face_vertex_num[idx_end - 1] - 1;
  vertex_num2 = face_vertex_num[idx_start] - 1;

  /* Get vertex coordinates */

  for (i_dim = 0; i_dim < dim; i_dim++)
    vect[i_dim] =  vertex_coords[vertex_num2*dim + i_dim]
                 - vertex_coords[vertex_num1*dim + i_dim];

  edge_cos_alpha = _COSINE_3D(vect, face_normal);
  edge_cos_alpha = CS_ABS(edge_cos_alpha);
  cos_alpha = CS_MIN(cos_alpha, edge_cos_alpha);

  *face_warping = 90. - acos(cos_alpha) * rad_to_deg;

}

/*----------------------------------------------------------------------------
 * Compute angle between face normal and segment dased on centers of the
 * adjacent cells. It evaluates a difference with orthogonality.
 *
 * parameters:
 *   mesh             --> pointer to mesh structure
 *   mesh_quantities  --> pointer to mesh quantities structures
 *   fac_warping      <-> array for internal faces
 *   fbr_warping      <-> array for border faces
 *----------------------------------------------------------------------------*/

static void
_compute_warping(const cs_maillage_t      *mesh,
                 const cs_maillage_grd_t  *mesh_quantities,
                 cs_real_t                 fac_warping[],
                 cs_real_t                 fbr_warping[])
{
  cs_int_t  i_fac, i_dim;
  cs_int_t  idx_start, idx_end;

  cs_real_t  face_normal[3];

  const cs_int_t  dim = mesh->dim;
  const cs_int_t  *fac_vtx_idx = mesh->pos_fac_som;
  const cs_int_t  *fbr_vtx_idx = mesh->pos_fbr_som;

  assert(dim == 3);

  /* Compute warping for internal faces */
  /*------------------------------------*/

  for (i_fac = 0; i_fac < mesh->nbr_fac; i_fac++) {

    /* Get normal to the face */

    for (i_dim = 0; i_dim < dim; i_dim++)
      face_normal[i_dim] = mesh_quantities->surf_fac[i_fac*dim + i_dim];

    /* Evaluate warping for each edge */

    idx_start = fac_vtx_idx[i_fac] - 1;
    idx_end = fac_vtx_idx[i_fac + 1] - 1;

    _get_face_warping(idx_start,
                      idx_end,
                      face_normal,
                      mesh->val_fac_som,
                      mesh->coo_som,
                      &(fac_warping[i_fac]));

  } /* End of loop on internal faces */

  /* Compute warping for border faces */
  /*----------------------------------*/

  for (i_fac = 0; i_fac < mesh->nbr_fbr; i_fac++) {

    /* Get face normal */

    for (i_dim = 0; i_dim < dim; i_dim++)
      face_normal[i_dim] = mesh_quantities->surf_fbr[i_fac*dim + i_dim];

    /* Evaluate warping for each edge */

    idx_start = fbr_vtx_idx[i_fac] - 1;
    idx_end = fbr_vtx_idx[i_fac + 1] - 1;

    _get_face_warping(idx_start,
                      idx_end,
                      face_normal,
                      mesh->val_fbr_som,
                      mesh->coo_som,
                      &(fbr_warping[i_fac]));

  } /* End of loop on border faces */

}

/*----------------------------------------------------------------------------
 * Transform face values to cell values using the maximum value
 * of a cell's faces.
 *
 * parameters:
 *   mesh          --> pointer to mesh structure
 *   default_value --> default value for initialization
 *   fac_val       --> interior face values
 *   fbr_val       --> boundary face values
 *   cel_val       <-- cell values
 *----------------------------------------------------------------------------*/

static void
_cell_from_max_face(const cs_maillage_t  *mesh,
                    const cs_real_t       default_value,
                    const cs_real_t       fac_val[],
                    const cs_real_t       fbr_val[],
                    cs_real_t             cel_val[])
{
  cs_int_t  ii, jj;
  cs_int_t  cell_id;

  /* Default initialization */

  for (ii = 0; ii < mesh->nbr_cel_etendu; ii++)
    cel_val[ii] = default_value;

  /* Distribution */

  if (fac_val != NULL) {

    for (ii = 0; ii < mesh->nbr_fac; ii++) {
      for (jj = 0; jj < 2; jj++) {
        cell_id = mesh->fac_cel[2*ii + jj] - 1;
        if (fac_val[ii] > cel_val[cell_id])
          cel_val[cell_id] = fac_val[ii];
      }
    }

  }

  if (fbr_val != NULL) {

    for (ii = 0; ii < mesh->nbr_fbr; ii++) {
      cell_id = mesh->fbr_cel[ii] - 1;
      if (fbr_val[ii] > cel_val[cell_id])
        cel_val[cell_id] = fbr_val[ii];
    }

  }
}

/*----------------------------------------------------------------------------
 * Transform face values to vertex values using the maximum value
 * of a vertices's connected faces.
 *
 * parameters:
 *   mesh          --> pointer to mesh structure
 *   default_value --> default value for initialization
 *   fac_val       --> interior face values
 *   fbr_val       --> boundary face values
 *   cel_val       <-- cell values
 *----------------------------------------------------------------------------*/

static void
_vtx_from_max_face(const cs_maillage_t  *mesh,
                   const cs_real_t       default_value,
                   const cs_real_t       fac_val[],
                   const cs_real_t       fbr_val[],
                   cs_real_t             vtx_val[])
{
  cs_int_t  ii, jj, idx_start, idx_end;
  cs_int_t  vtx_id;

  /* Default initialization */

  for (ii = 0; ii < mesh->nbr_som; ii++)
    vtx_val[ii] = default_value;

  /* Distribution */

  if (fac_val != NULL && mesh->pos_fac_som != NULL) {

    for (ii = 0; ii < mesh->nbr_fac; ii++) {

      idx_start = mesh->pos_fac_som[ii] - 1;
      idx_end = mesh->pos_fac_som[ii + 1] - 1;

      for (jj = idx_start; jj < idx_end; jj++) {
        vtx_id = mesh->val_fac_som[jj] - 1;
        if (fac_val[ii] > vtx_val[vtx_id])
          vtx_val[vtx_id] = fac_val[ii];
      }
    }

  }

  if (fbr_val != NULL && mesh->pos_fbr_som != NULL) {

    for (ii = 0; ii < mesh->nbr_fbr; ii++) {

      idx_start = mesh->pos_fbr_som[ii] - 1;
      idx_end = mesh->pos_fbr_som[ii + 1] - 1;

      for (jj = idx_start; jj < idx_end; jj++) {
        vtx_id = mesh->val_fbr_som[jj] - 1;
        if (fbr_val[ii] > vtx_val[vtx_id])
          vtx_val[vtx_id] = fbr_val[ii];
      }
    }

  }
}

/*============================================================================
 * Public function definitions for Fortran API
 *============================================================================*/

/*============================================================================
 * Public function definitions
 *============================================================================*/

/*----------------------------------------------------------------------------
 * Compute mesh quality indicators
 *
 * parameters:
 *   mesh             --> pointer to mesh structure.
 *   mesh_quantities  --> pointer to mesh quantities structure.
 *----------------------------------------------------------------------------*/

void
cs_mesh_quality(const cs_maillage_t      *mesh,
                const cs_maillage_grd_t  *mesh_quantities)
{
  cs_int_t  i;

  cs_bool_t compute_volume = CS_TRUE;
  cs_bool_t compute_weighting = CS_TRUE;
  cs_bool_t compute_orthogonality = CS_TRUE;
  cs_bool_t compute_warping = CS_TRUE;

  cs_bool_t vol_fields = CS_FALSE;
  cs_bool_t brd_fields = CS_FALSE;

  cs_real_t *face_to_cell = NULL;
  cs_real_t *face_to_vtx = NULL;

  double  *working_array = NULL;

  const cs_int_t  n_vertices = mesh->nbr_som;
  const cs_int_t  n_fac = mesh->nbr_fac;
  const cs_int_t  n_fbr = mesh->nbr_fbr;
  const cs_int_t  n_cells = mesh->nbr_cel;
  const cs_int_t  n_cells_ext = mesh->nbr_cel_etendu;

  /* Check input data */

  assert(mesh_quantities->surf_fac != NULL);
  assert(mesh_quantities->cdg_fac != NULL);
  assert(mesh_quantities->cen_cel != NULL);
  assert(mesh_quantities->vol_cel != NULL);

  /* Determine if resulting fields should be exported on the volume
     and border meshes (depending on their existence); */

  /* Note that n_vertices or n_cells should never be zero on any
     rank (unlike n_fbr), so if face_to_cell/face_to_vtx is allocated
     on any rank, it should be allocated on all ranks;
     We can thus use this pointer for tests safely */

  /* TODO:
     define an option to distribute face values to cells, vertices, or both */

  if (cs_post_existe_maillage(-1)) {
    vol_fields = CS_TRUE;
    BFT_MALLOC(face_to_cell, CS_MAX(n_cells_ext, n_vertices), cs_real_t);
    face_to_vtx = face_to_cell;
  }

  if (cs_post_existe_maillage(-2))
    brd_fields = CS_TRUE;

  /* TODO
     For the moment, we export the mesh at this stage; this should be moved
     once PSTEMA has been moved from CALTRI to an earlier step. */

  cs_post_ecrit_maillages(1,0);

  cs_post_activer_writer(0, 1);

  /* Evaluate mesh quality criteria */
  /*--------------------------------*/

  /*--------------*/
  /* Face warping */
  /*--------------*/

  if (compute_warping == CS_TRUE) {

    double  *fac_warping = NULL, *fbr_warping = NULL;

    BFT_MALLOC(working_array, n_fac + n_fbr, double);

    for (i = 0; i < n_fac + n_fbr; i++)
      working_array[i] = 0.;

    fac_warping = working_array;
    fbr_warping = working_array + n_fac;

    _compute_warping(mesh,
                     mesh_quantities,
                     fac_warping,
                     fbr_warping);

    /* Display histograms */

    bft_printf(_("\n  Histogramme du gauchissement des faces internes :\n\n"));

    _display_histograms(n_fac, fac_warping);

    if (mesh->nbr_fbr_glob > 0) {

      bft_printf(_("\n  Histogramme du gauchissement des faces de bord:\n\n"));

      _display_histograms(n_fbr, fbr_warping);

    }

    /* Post processing */

    if (vol_fields == CS_TRUE) {

      if (face_to_cell != NULL) {
        _cell_from_max_face(mesh, 0., fac_warping, fbr_warping, face_to_cell);
        cs_post_ecrit_var(-1,
                          "Face_Warp_c_max",
                          1,
                          CS_FALSE,
                          CS_TRUE,
                          CS_POST_TYPE_cs_real_t,
                          -1,
                          0.0,
                          face_to_cell,
                          NULL,
                          NULL);
      }
      if (face_to_vtx != NULL) {
        _vtx_from_max_face(mesh, 0., fac_warping, fbr_warping, face_to_vtx);
        cs_post_ecrit_var_som(-1,
                              "Face_Warp_v_max",
                              1,
                              CS_FALSE,
                              CS_TRUE,
                              CS_POST_TYPE_cs_real_t,
                              -1,
                              0.0,
                              face_to_vtx);
      }

    }

    if (brd_fields == CS_TRUE)
        cs_post_ecrit_var(-2,
                          "Face_Warp",
                          1,
                          CS_FALSE,
                          CS_TRUE,
                          CS_POST_TYPE_cs_real_t,
                          -1,
                          0.0,
                          NULL,
                          NULL,
                          fbr_warping);

    BFT_FREE(working_array);

  }

  /*----------------------------------------------*/
  /* Weighting and center offsetting coefficients */
  /*----------------------------------------------*/

  if (compute_weighting == CS_TRUE) {

    double  *weighting = NULL, *offsetting = NULL;

    /* Only defined on internal faces */

    BFT_MALLOC(working_array, 2*n_fac, double);

    for (i = 0; i < 2*n_fac; i++)
      working_array[i] = 0.;

    weighting = working_array;
    offsetting = working_array + n_fac;

    _compute_weighting_offsetting(mesh,
                                  mesh_quantities,
                                  weighting,
                                  offsetting);

    /* Display histograms */

    bft_printf(_("\n  Histogramme des coefficients de pondération des"
               " faces internes :\n\n"));

    _display_histograms(n_fac, weighting);

    bft_printf(_("\n  Histogramme des coefficients de décentrement"
               " des faces internes :\n\n"));

    _display_histograms(n_fac, offsetting);

    /* Post processing */

    if (vol_fields == CS_TRUE) {

      if (face_to_cell != NULL) {
        _cell_from_max_face(mesh, 0.5, weighting, NULL, face_to_cell);
        cs_post_ecrit_var(-1,
                          "Weighting_c_max",
                          1,
                          CS_FALSE,
                          CS_TRUE,
                          CS_POST_TYPE_cs_real_t,
                          -1,
                          0.0,
                          face_to_cell,
                          NULL,
                          NULL);
      }
      if (face_to_vtx != NULL) {
        _vtx_from_max_face(mesh, 0.5, weighting, NULL, face_to_vtx);
        cs_post_ecrit_var_som(-1,
                              "Weighting_v_max",
                              1,
                              CS_FALSE,
                              CS_TRUE,
                              CS_POST_TYPE_cs_real_t,
                              -1,
                              0.0,
                              face_to_vtx);
      }

      if (face_to_cell != NULL) {
        _cell_from_max_face(mesh, 0., offsetting, NULL, face_to_cell);
        cs_post_ecrit_var(-1,
                          "Offset_c_max",
                          1,
                          CS_FALSE,
                          CS_TRUE,
                          CS_POST_TYPE_cs_real_t,
                          -1,
                          0.0,
                          face_to_cell,
                          NULL,
                          NULL);
      }
      if (face_to_vtx != NULL) {
        _vtx_from_max_face(mesh, 0., offsetting, NULL, face_to_vtx);
        cs_post_ecrit_var_som(-1,
                              "Offset_v_max",
                              1,
                              CS_FALSE,
                              CS_TRUE,
                              CS_POST_TYPE_cs_real_t,
                              -1,
                              0.0,
                              face_to_vtx);
      }

    }

    BFT_FREE(working_array);

  }

  /*---------------------*/
  /* Angle orthogonality */
  /*---------------------*/

  if (compute_orthogonality == CS_TRUE) {

    double  *fac_ortho = NULL, *fbr_ortho = NULL;

    BFT_MALLOC(working_array, n_fac + n_fbr, double);

    for (i = 0; i < n_fac + n_fbr; i++)
      working_array[i] = 0.;

    fac_ortho = working_array;
    fbr_ortho = working_array + n_fac;

    _compute_orthogonality(mesh,
                           mesh_quantities,
                           fac_ortho,
                           fbr_ortho);

    /* Display histograms */

    bft_printf(_("\n  Histogramme de la non-orthogonalité des"
               " faces internes (en degrés) :\n\n"));

    _display_histograms(n_fac, fac_ortho);

    if (mesh->nbr_fbr_glob > 0) {

      bft_printf(_("\n  Histogramme de la non-orthogonalité des"
                 " faces de bord (en degrés) :\n\n"));

      _display_histograms(n_fbr, fbr_ortho);

    }

    /* Post processing */

    if (vol_fields == CS_TRUE) {

      if (face_to_cell != NULL) {
        _cell_from_max_face(mesh, 0., fac_ortho, fbr_ortho, face_to_cell);
        cs_post_ecrit_var(-1,
                          "Non_Ortho_c_max",
                          1,
                          CS_FALSE,
                          CS_TRUE,
                          CS_POST_TYPE_cs_real_t,
                          -1,
                          0.0,
                          face_to_cell,
                          NULL,
                          NULL);
      }
      if (face_to_vtx != NULL) {
        _vtx_from_max_face(mesh, 0., fac_ortho, fbr_ortho, face_to_vtx);
        cs_post_ecrit_var_som(-1,
                              "Non_Ortho_v_max",
                              1,
                              CS_FALSE,
                              CS_TRUE,
                              CS_POST_TYPE_cs_real_t,
                              -1,
                              0.0,
                              face_to_vtx);
      }

    }

    if (brd_fields == CS_TRUE)
        cs_post_ecrit_var(-2,
                          "Non_Ortho",
                          1,
                          CS_FALSE,
                          CS_TRUE,
                          CS_POST_TYPE_cs_real_t,
                          -1,
                          0.0,
                          NULL,
                          NULL,
                          fbr_ortho);

    BFT_FREE(working_array);

  }

  if (face_to_cell != NULL)
    BFT_FREE(face_to_cell);

  /*-------------*/
  /* Cell volume */
  /*-------------*/

  if (compute_volume == CS_TRUE) {

    /* Display histograms */

    bft_printf(_("\n  Histogramme du volume des cellules :\n\n")) ;

    _display_histograms(n_cells, mesh_quantities->vol_cel);

    /* Post processing */

    if (vol_fields == CS_TRUE)
        cs_post_ecrit_var(-1,
                          "Cell_Volume",
                          1,
                          CS_FALSE,
                          CS_TRUE,
                          CS_POST_TYPE_cs_real_t,
                          -1,
                          0.0,
                          mesh_quantities->vol_cel,
                          NULL,
                          NULL);

  }

}


/*----------------------------------------------------------------------------*/

/* Delete local macros */

#undef _CROSS_PRODUCT_3D
#undef _DOT_PRODUCT_3D
#undef _MODULE_3D
#undef _COSINE_3D

#ifdef __cplusplus
}
#endif /* __cplusplus */
