/*============================================================================
*
* 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 */
syntax highlighted by Code2HTML, v. 0.9.1