1
/*============================================================================
3
* This file is part of the Code_Saturne Kernel, element of the
4
* Code_Saturne CFD tool.
6
* Copyright (C) 1998-2009 EDF S.A., France
8
* contact: saturne-support@edf.fr
10
* The Code_Saturne Kernel is free software; you can redistribute it
11
* and/or modify it under the terms of the GNU General Public License
12
* as published by the Free Software Foundation; either version 2 of
13
* the License, or (at your option) any later version.
15
* The Code_Saturne Kernel is distributed in the hope that it will be
16
* useful, but WITHOUT ANY WARRANTY; without even the implied warranty
17
* of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18
* GNU General Public License for more details.
20
* You should have received a copy of the GNU General Public License
21
* along with the Code_Saturne Kernel; if not, write to the
22
* Free Software Foundation, Inc.,
23
* 51 Franklin St, Fifth Floor,
24
* Boston, MA 02110-1301 USA
26
*============================================================================*/
28
/*============================================================================
29
* Functions checking the coherency of the mesh
30
*============================================================================*/
32
#if defined(HAVE_CONFIG_H)
33
#include "cs_config.h"
36
/*----------------------------------------------------------------------------
37
* Standard C library headers
38
*----------------------------------------------------------------------------*/
44
/*----------------------------------------------------------------------------
46
*----------------------------------------------------------------------------*/
49
#include <bft_error.h>
50
#include <bft_printf.h>
52
/*----------------------------------------------------------------------------
54
*----------------------------------------------------------------------------*/
56
/*----------------------------------------------------------------------------
58
*----------------------------------------------------------------------------*/
62
#include "cs_mesh_quantities.h"
65
/*----------------------------------------------------------------------------
66
* Header for the current file
67
*----------------------------------------------------------------------------*/
69
#include "cs_mesh_coherency.h"
71
/*----------------------------------------------------------------------------*/
75
/*============================================================================
76
* Local structure and type definitions
77
*============================================================================*/
79
/*=============================================================================
80
* Local Macro definitions
81
*============================================================================*/
83
#define CS_MESH_COHERENCY_TOLERANCE 0.05
85
/*============================================================================
86
* Global static variables
87
*============================================================================*/
89
/*=============================================================================
90
* Private function definitions
91
*============================================================================*/
93
/*----------------------------------------------------------------------------
94
* Check the coherency of the internal face -> cells connectivity.
95
*----------------------------------------------------------------------------*/
102
const cs_mesh_t *mesh = cs_glob_mesh;
104
bft_printf(_(" Checking the face -> cells connectivity coherency\n"));
106
for (i = 0; i < mesh->n_i_faces; i++) {
108
if (mesh->i_face_cells[2*i] == 0 || mesh->i_face_cells[2*i+1] == 0)
109
bft_error(__FILE__, __LINE__, 0,
110
_("Internal face -> cells connectivity value not initialized\n"
111
"for face: %d (cell_num1 = %d and cell_num2 = %d)\n"),
112
i+1, mesh->i_face_cells[2*i], mesh->i_face_cells[2*i+1]);
118
/*=============================================================================
119
* Public function definitions
120
*============================================================================*/
122
/*----------------------------------------------------------------------------
123
* Check the coherency of the global mesh structure.
124
*----------------------------------------------------------------------------*/
127
cs_mesh_coherency_check(void)
129
cs_int_t i, j, coord_id, cell_id, face_id, vtx_id;
130
cs_real_t test, delta_mean, delta_neighbor;
131
cs_real_t _min, _max, coord;
133
cs_real_t delta_mean_mult = 1.0;
135
cs_real_t *minmax_buffer = NULL, *compute_buffer = NULL;
136
cs_real_t *min = NULL, *max = NULL, *mean = NULL, *delta = NULL;
138
const cs_mesh_t *mesh = cs_glob_mesh;
139
const cs_mesh_quantities_t *mesh_quantities = cs_glob_mesh_quantities;
140
const cs_int_t n_cells = mesh->n_cells;
141
const cs_int_t n_cells_with_ghosts = mesh->n_cells_with_ghosts;
143
const cs_int_t *ifacel = mesh->i_face_cells;
144
const cs_int_t *bfacel = mesh->b_face_cells;
145
const cs_real_t *vtx_coord = mesh->vtx_coord;
147
bft_printf(_("\n Checking the mesh structure coherency:\n"));
149
/* Check internal face -> cells connectivity coherency */
153
/* allocate and initialize buffers */
155
BFT_MALLOC(minmax_buffer, 2*n_cells, cs_real_t);
156
BFT_MALLOC(compute_buffer, 6*n_cells_with_ghosts, cs_real_t);
159
max = minmax_buffer + n_cells;
160
mean = compute_buffer;
161
delta = compute_buffer + 3*n_cells_with_ghosts;
163
/* Loop on coordinates */
165
bft_printf(_(" Coherency criteria definition\n"));
167
for (coord_id = 0; coord_id < 3; coord_id++) {
169
for (i = 0; i < n_cells; i++) {
174
for (i = 0; i < n_cells_with_ghosts; i++) {
175
mean[3*i + coord_id] = DBL_MIN;
176
delta[3*i + coord_id] = DBL_MIN;
179
/* Loop on internal faces */
181
for (face_id = 0; face_id < mesh->n_i_faces; face_id++) {
183
const cs_int_t *face_vtx_idx = mesh->i_face_vtx_idx;
184
const cs_int_t *face_vtx_lst = mesh->i_face_vtx_lst;
189
for (i = face_vtx_idx[face_id]-1; i < face_vtx_idx[face_id+1]-1; i++) {
191
vtx_id = face_vtx_lst[i]-1;
192
coord = vtx_coord[3*vtx_id + coord_id];
193
_min = CS_MIN(_min, coord);
194
_max = CS_MAX(_max, coord);
198
for (j = 0; j < 2; j++) {
200
cell_id = ifacel[2*face_id + j] - 1;
202
if (cell_id < n_cells) {
203
max[cell_id] = CS_MAX(max[cell_id], _max);
204
min[cell_id] = CS_MIN(min[cell_id], _min);
209
} /* End of loop on internal faces */
211
/* Loop on border faces */
213
for (face_id = 0; face_id < mesh->n_b_faces; face_id++) {
215
const cs_int_t *face_vtx_idx = mesh->b_face_vtx_idx;
216
const cs_int_t *face_vtx_lst = mesh->b_face_vtx_lst;
221
for (i = face_vtx_idx[face_id]-1; i < face_vtx_idx[face_id+1]-1; i++) {
223
vtx_id = face_vtx_lst[i]-1;
224
coord = vtx_coord[3*vtx_id + coord_id];
225
_min = CS_MIN(_min, coord);
226
_max = CS_MAX(_max, coord);
230
cell_id = bfacel[face_id] - 1;
232
assert(cell_id < n_cells);
233
max[cell_id] = CS_MAX(max[cell_id], _max);
234
min[cell_id] = CS_MIN(min[cell_id], _min);
236
} /* End of loop on border faces */
238
/* Compute mean and delta for this coordinate */
240
for (cell_id = 0; cell_id < n_cells; cell_id++) {
242
mean[3*cell_id + coord_id] = (max[cell_id] + min[cell_id])/2;
243
delta[3*cell_id + coord_id] = max[cell_id] - min[cell_id];
245
assert(delta[3*cell_id + coord_id] > 0);
249
} /* End of loop on coordinates */
251
/* Synchronize variable */
253
if (mesh->halo != NULL) {
255
cs_halo_sync_var_strided(mesh->halo, mesh->halo_type, mean, 3);
256
cs_halo_sync_var_strided(mesh->halo, mesh->halo_type, delta, 3);
260
/* Synchronization for periodicity */
262
if (mesh->n_init_perio > 0) {
264
cs_real_t *delta_buffer;
266
BFT_MALLOC(delta_buffer, 3*n_cells_with_ghosts, cs_real_t);
268
/* De-interlace delta arrays for periodicity exchange;
269
Also add factor of 1.8 (> sqrt(3)) as the longest diagonal
270
of a box of side 1 is sqrt(3), and may find itself aligned
271
with axes after rotation. */
273
if (mesh->have_rotation_perio == 1) {
275
delta_mean_mult = 1.0/1.8;
277
for (cell_id = 0; cell_id < n_cells_with_ghosts; cell_id++) {
278
cs_real_t delta_max = delta[3*cell_id];
279
if (delta[3*cell_id + 1] > delta_max)
280
delta_max = delta[3*cell_id + 1];
281
if (delta[3*cell_id + 2] > delta_max)
282
delta_max = delta[3*cell_id + 2];
284
delta_buffer[ cell_id] = delta_max;
285
delta_buffer[ n_cells_with_ghosts + cell_id] = delta_max;
286
delta_buffer[2*n_cells_with_ghosts + cell_id] = delta_max;
292
for (coord_id = 0; coord_id < 3; coord_id++) {
293
for (cell_id = 0; cell_id < n_cells_with_ghosts; cell_id++)
294
delta_buffer[coord_id*n_cells_with_ghosts + cell_id]
295
= delta[3*cell_id + coord_id];
300
cs_perio_sync_var_vect(mesh->halo,
304
delta_buffer + n_cells_with_ghosts,
305
delta_buffer + 2*n_cells_with_ghosts);
307
for (coord_id = 0; coord_id < 3; coord_id++) {
308
for (cell_id = n_cells; cell_id < n_cells_with_ghosts; cell_id++)
309
delta[3*cell_id + coord_id] =
310
delta_buffer[coord_id*n_cells_with_ghosts + cell_id];
313
BFT_FREE(delta_buffer);
315
cs_perio_sync_coords(mesh->halo, mesh->halo_type, mean);
319
for (coord_id = 0; coord_id < 3; coord_id++) {
321
bft_printf(_(" Coherency verification on coordinates %d\n"),
324
/* Test coherency on the standard neighborhood */
326
for (face_id = 0; face_id < mesh->n_i_faces; face_id++) {
328
cs_int_t cell_id1 = ifacel[2*face_id] - 1;
329
cs_int_t cell_id2 = ifacel[2*face_id + 1] - 1;
330
cs_real_t delta1 = CS_ABS(delta[3*cell_id1 + coord_id]);
331
cs_real_t delta2 = CS_ABS(delta[3*cell_id2 + coord_id]);
332
cs_real_t mean1 = mean[3*cell_id1 + coord_id];
333
cs_real_t mean2 = mean[3*cell_id2 + coord_id];
335
delta_mean = CS_ABS(mean2 - mean1)*delta_mean_mult;
336
delta_neighbor = (delta1 + delta2)/2.;
338
test = (1 + CS_MESH_COHERENCY_TOLERANCE)*delta_neighbor - delta_mean;
342
cs_real_t *cell_cen = mesh_quantities->cell_cen;
344
bft_printf(_("\nInfo on cell 1: %d\n"
345
" cell center: %12.3g %12.3g %12.3g\n"
347
" box center: %12.3g\n"),
348
cell_id1+1, cell_cen[3*cell_id1], cell_cen[3*cell_id1+1],
349
cell_cen[3*cell_id1+2], delta1, mean1);
351
bft_printf(_("\nInfo on cell 2: %d\n"
352
" cell center: %12.3g %12.3g %12.3g\n"
354
" box center: %12.3g\n"),
355
cell_id2+1, cell_cen[3*cell_id2], cell_cen[3*cell_id2+1],
356
cell_cen[3*cell_id2+2], delta2, mean2);
359
bft_error(__FILE__, __LINE__, 0,
360
_("\nCoherency error in standard halo\n"
361
"between cells %d and %d: test = %g\n"
362
"(delta = %g, delta_mean = %g)\n"),
363
cell_id1+1, cell_id2+1, test, delta_neighbor, delta_mean);
367
} /* End of loop on internal faces */
369
if (mesh->cell_cells_idx != NULL) {
371
cs_int_t *cell_cells_idx = mesh->cell_cells_idx;
373
for (cell_id = 0; cell_id < n_cells; cell_id++) {
375
for (i = cell_cells_idx[cell_id]-1;
376
i < cell_cells_idx[cell_id+1]-1; i++) {
378
cs_int_t cell_id2 = mesh->cell_cells_lst[i] - 1;
379
cs_real_t delta1 = CS_ABS(delta[3*cell_id + coord_id]);
380
cs_real_t delta2 = CS_ABS(delta[3*cell_id2 + coord_id]);
381
cs_real_t mean1 = mean[3*cell_id + coord_id];
382
cs_real_t mean2 = mean[3*cell_id2 + coord_id];
384
delta_mean = CS_ABS(mean2 - mean1)*delta_mean_mult;
385
delta_neighbor = (delta1 + delta2)/2.;
387
test = (1 + CS_MESH_COHERENCY_TOLERANCE)*delta_neighbor - delta_mean;
391
cs_real_t *cell_cen = mesh_quantities->cell_cen;
393
bft_printf(_("\nInfo on cell 1: %d\n"
394
" cell center: %12.3g %12.3g %12.3g\n"
396
" box center: %12.3g\n"),
397
cell_id+1, cell_cen[3*cell_id], cell_cen[3*cell_id+1],
398
cell_cen[3*cell_id+2], delta1, mean1);
400
bft_printf(_("\nInfo on cell 2: %d\n"
401
" cell center: %12.3g %12.3g %12.3g\n"
403
" box center: %12.3g\n"),
404
cell_id2+1, cell_cen[3*cell_id2], cell_cen[3*cell_id2+1],
405
cell_cen[3*cell_id2+2], delta2, mean2);
408
bft_error(__FILE__, __LINE__, 0,
409
_("\nCoherency error in extended halo\n"
410
"between cells %d and %d: test = %g\n"
411
"(delta = %g, delta_mean = %g)\n"),
412
cell_id+1, cell_id2+1, test, delta_neighbor, delta_mean);
418
} /* End of loop on cells */
420
} /* End of treatment of the exetended neighborhood */
422
} /* End of loop on coordinates */
426
BFT_FREE(compute_buffer);
427
BFT_FREE(minmax_buffer);
429
bft_printf(_(" End of coherency check of the mesh structure.\n"));
433
/*----------------------------------------------------------------------------*/