/*============================================================================ * * 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 #include #include #include #include /*---------------------------------------------------------------------------- * BFT library headers *----------------------------------------------------------------------------*/ #include #include #include /*---------------------------------------------------------------------------- * 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 */