~ubuntu-branches/ubuntu/precise/code-saturne/precise

« back to all changes in this revision

Viewing changes to src/base/cs_mesh_coherency.c

  • Committer: Package Import Robot
  • Author(s): Sylvestre Ledru
  • Date: 2011-11-01 17:43:32 UTC
  • mto: (6.1.7 sid)
  • mto: This revision was merged to the branch mainline in revision 11.
  • Revision ID: package-import@ubuntu.com-20111101174332-tl4vk45no0x3emc3
Tags: upstream-2.1.0
ImportĀ upstreamĀ versionĀ 2.1.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
/*============================================================================
2
 
 *
3
 
 *     This file is part of the Code_Saturne Kernel, element of the
4
 
 *     Code_Saturne CFD tool.
5
 
 *
6
 
 *     Copyright (C) 1998-2009 EDF S.A., France
7
 
 *
8
 
 *     contact: saturne-support@edf.fr
9
 
 *
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.
14
 
 *
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.
19
 
 *
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
25
 
 *
26
 
 *============================================================================*/
27
 
 
28
 
/*============================================================================
29
 
 * Functions checking the coherency of the mesh
30
 
 *============================================================================*/
31
 
 
32
 
#if defined(HAVE_CONFIG_H)
33
 
#include "cs_config.h"
34
 
#endif
35
 
 
36
 
/*----------------------------------------------------------------------------
37
 
 * Standard C library headers
38
 
 *----------------------------------------------------------------------------*/
39
 
 
40
 
#include <stdio.h>
41
 
#include <float.h>
42
 
#include <assert.h>
43
 
 
44
 
/*----------------------------------------------------------------------------
45
 
 * BFT library headers
46
 
 *----------------------------------------------------------------------------*/
47
 
 
48
 
#include <bft_mem.h>
49
 
#include <bft_error.h>
50
 
#include <bft_printf.h>
51
 
 
52
 
/*----------------------------------------------------------------------------
53
 
 * FVM library headers
54
 
 *----------------------------------------------------------------------------*/
55
 
 
56
 
/*----------------------------------------------------------------------------
57
 
 *  Local headers
58
 
 *----------------------------------------------------------------------------*/
59
 
 
60
 
#include "cs_halo.h"
61
 
#include "cs_mesh.h"
62
 
#include "cs_mesh_quantities.h"
63
 
#include "cs_perio.h"
64
 
 
65
 
/*----------------------------------------------------------------------------
66
 
 *  Header for the current file
67
 
 *----------------------------------------------------------------------------*/
68
 
 
69
 
#include "cs_mesh_coherency.h"
70
 
 
71
 
/*----------------------------------------------------------------------------*/
72
 
 
73
 
BEGIN_C_DECLS
74
 
 
75
 
/*============================================================================
76
 
 * Local structure and type definitions
77
 
 *============================================================================*/
78
 
 
79
 
/*=============================================================================
80
 
 * Local Macro definitions
81
 
 *============================================================================*/
82
 
 
83
 
#define CS_MESH_COHERENCY_TOLERANCE  0.05
84
 
 
85
 
/*============================================================================
86
 
 *  Global static variables
87
 
 *============================================================================*/
88
 
 
89
 
/*=============================================================================
90
 
 * Private function definitions
91
 
 *============================================================================*/
92
 
 
93
 
/*----------------------------------------------------------------------------
94
 
 * Check the coherency of the internal face -> cells connectivity.
95
 
 *----------------------------------------------------------------------------*/
96
 
 
97
 
static void
98
 
_check_ifacel(void)
99
 
{
100
 
  cs_int_t  i;
101
 
 
102
 
  const cs_mesh_t  *mesh = cs_glob_mesh;
103
 
 
104
 
  bft_printf(_("    Checking the face -> cells connectivity coherency\n"));
105
 
 
106
 
  for (i = 0; i < mesh->n_i_faces; i++) {
107
 
 
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]);
113
 
 
114
 
  }
115
 
 
116
 
}
117
 
 
118
 
/*=============================================================================
119
 
 * Public function definitions
120
 
 *============================================================================*/
121
 
 
122
 
/*----------------------------------------------------------------------------
123
 
 * Check the coherency of the global mesh structure.
124
 
 *----------------------------------------------------------------------------*/
125
 
 
126
 
void
127
 
cs_mesh_coherency_check(void)
128
 
{
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;
132
 
 
133
 
  cs_real_t delta_mean_mult = 1.0;
134
 
 
135
 
  cs_real_t  *minmax_buffer = NULL, *compute_buffer = NULL;
136
 
  cs_real_t  *min = NULL, *max = NULL, *mean = NULL, *delta = NULL;
137
 
 
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;
142
 
 
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;
146
 
 
147
 
  bft_printf(_("\n Checking the mesh structure coherency:\n"));
148
 
 
149
 
  /* Check internal face -> cells connectivity coherency */
150
 
 
151
 
  _check_ifacel();
152
 
 
153
 
  /* allocate and initialize buffers */
154
 
 
155
 
  BFT_MALLOC(minmax_buffer, 2*n_cells, cs_real_t);
156
 
  BFT_MALLOC(compute_buffer, 6*n_cells_with_ghosts, cs_real_t);
157
 
 
158
 
  min = minmax_buffer;
159
 
  max = minmax_buffer + n_cells;
160
 
  mean = compute_buffer;
161
 
  delta = compute_buffer + 3*n_cells_with_ghosts;
162
 
 
163
 
  /* Loop on coordinates */
164
 
 
165
 
  bft_printf(_("    Coherency criteria definition\n"));
166
 
 
167
 
  for (coord_id = 0; coord_id < 3; coord_id++) {
168
 
 
169
 
    for (i = 0; i < n_cells; i++) {
170
 
      min[i] = DBL_MAX;
171
 
      max[i] =-DBL_MAX;
172
 
    }
173
 
 
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;
177
 
    }
178
 
 
179
 
    /* Loop on internal faces */
180
 
 
181
 
    for (face_id = 0; face_id < mesh->n_i_faces; face_id++) {
182
 
 
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;
185
 
 
186
 
      _min = DBL_MAX;
187
 
      _max =-DBL_MAX;
188
 
 
189
 
      for (i = face_vtx_idx[face_id]-1; i < face_vtx_idx[face_id+1]-1; i++) {
190
 
 
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);
195
 
 
196
 
      }
197
 
 
198
 
      for (j = 0; j < 2; j++) {
199
 
 
200
 
        cell_id = ifacel[2*face_id + j] - 1;
201
 
 
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);
205
 
        }
206
 
 
207
 
      }
208
 
 
209
 
    } /* End of loop on internal faces */
210
 
 
211
 
    /* Loop on border faces */
212
 
 
213
 
    for (face_id = 0; face_id < mesh->n_b_faces; face_id++) {
214
 
 
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;
217
 
 
218
 
      _min = DBL_MAX;
219
 
      _max =-DBL_MAX;
220
 
 
221
 
      for (i = face_vtx_idx[face_id]-1; i < face_vtx_idx[face_id+1]-1; i++) {
222
 
 
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);
227
 
 
228
 
      }
229
 
 
230
 
      cell_id = bfacel[face_id] - 1;
231
 
 
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);
235
 
 
236
 
    } /* End of loop on border faces */
237
 
 
238
 
    /* Compute mean and delta for this coordinate */
239
 
 
240
 
    for (cell_id = 0; cell_id < n_cells; cell_id++) {
241
 
 
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];
244
 
 
245
 
      assert(delta[3*cell_id + coord_id] > 0);
246
 
 
247
 
    }
248
 
 
249
 
  } /* End of loop on coordinates */
250
 
 
251
 
  /* Synchronize variable */
252
 
 
253
 
  if (mesh->halo != NULL) {
254
 
 
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);
257
 
 
258
 
  }
259
 
 
260
 
  /* Synchronization for periodicity */
261
 
 
262
 
  if (mesh->n_init_perio > 0) {
263
 
 
264
 
    cs_real_t *delta_buffer;
265
 
 
266
 
    BFT_MALLOC(delta_buffer, 3*n_cells_with_ghosts, cs_real_t);
267
 
 
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. */
272
 
 
273
 
    if (mesh->have_rotation_perio == 1) {
274
 
 
275
 
      delta_mean_mult = 1.0/1.8;
276
 
 
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];
283
 
        delta_max *= 1.8;
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;
287
 
      }
288
 
 
289
 
    }
290
 
    else {
291
 
 
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];
296
 
      }
297
 
 
298
 
    }
299
 
 
300
 
    cs_perio_sync_var_vect(mesh->halo,
301
 
                           mesh->halo_type,
302
 
                           CS_PERIO_ROTA_COPY,
303
 
                           delta_buffer,
304
 
                           delta_buffer +   n_cells_with_ghosts,
305
 
                           delta_buffer + 2*n_cells_with_ghosts);
306
 
 
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];
311
 
    }
312
 
 
313
 
    BFT_FREE(delta_buffer);
314
 
 
315
 
    cs_perio_sync_coords(mesh->halo, mesh->halo_type, mean);
316
 
 
317
 
  }
318
 
 
319
 
  for (coord_id = 0; coord_id < 3; coord_id++) {
320
 
 
321
 
    bft_printf(_("    Coherency verification on coordinates %d\n"),
322
 
               coord_id+1);
323
 
 
324
 
    /* Test coherency on the standard neighborhood */
325
 
 
326
 
    for (face_id = 0; face_id < mesh->n_i_faces; face_id++) {
327
 
 
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];
334
 
 
335
 
      delta_mean = CS_ABS(mean2 - mean1)*delta_mean_mult;
336
 
      delta_neighbor = (delta1 + delta2)/2.;
337
 
 
338
 
      test = (1 + CS_MESH_COHERENCY_TOLERANCE)*delta_neighbor - delta_mean;
339
 
 
340
 
      if (test < 0) {
341
 
 
342
 
        cs_real_t  *cell_cen = mesh_quantities->cell_cen;
343
 
 
344
 
        bft_printf(_("\nInfo on cell 1: %d\n"
345
 
                     " cell center: %12.3g %12.3g %12.3g\n"
346
 
                     " delta:       %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);
350
 
 
351
 
        bft_printf(_("\nInfo on cell 2: %d\n"
352
 
                     " cell center: %12.3g %12.3g %12.3g\n"
353
 
                     " delta:       %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);
357
 
        bft_printf_flush();
358
 
 
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);
364
 
 
365
 
      }
366
 
 
367
 
    } /* End of loop on internal faces */
368
 
 
369
 
    if (mesh->cell_cells_idx != NULL) {
370
 
 
371
 
      cs_int_t  *cell_cells_idx = mesh->cell_cells_idx;
372
 
 
373
 
      for (cell_id = 0; cell_id < n_cells; cell_id++) {
374
 
 
375
 
        for (i = cell_cells_idx[cell_id]-1;
376
 
             i < cell_cells_idx[cell_id+1]-1; i++) {
377
 
 
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];
383
 
 
384
 
          delta_mean = CS_ABS(mean2 - mean1)*delta_mean_mult;
385
 
          delta_neighbor = (delta1 + delta2)/2.;
386
 
 
387
 
          test = (1 + CS_MESH_COHERENCY_TOLERANCE)*delta_neighbor - delta_mean;
388
 
 
389
 
          if (test < 0) {
390
 
 
391
 
            cs_real_t  *cell_cen = mesh_quantities->cell_cen;
392
 
 
393
 
            bft_printf(_("\nInfo on cell 1: %d\n"
394
 
                         " cell center: %12.3g %12.3g %12.3g\n"
395
 
                         " delta:       %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);
399
 
 
400
 
            bft_printf(_("\nInfo on cell 2: %d\n"
401
 
                         " cell center: %12.3g %12.3g %12.3g\n"
402
 
                         " delta:       %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);
406
 
            bft_printf_flush();
407
 
 
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);
413
 
 
414
 
          }
415
 
 
416
 
        }
417
 
 
418
 
      } /* End of loop on cells */
419
 
 
420
 
    } /* End of treatment of the exetended neighborhood */
421
 
 
422
 
  } /* End of loop on coordinates */
423
 
 
424
 
  /* Free memory */
425
 
 
426
 
  BFT_FREE(compute_buffer);
427
 
  BFT_FREE(minmax_buffer);
428
 
 
429
 
  bft_printf(_(" End of coherency check of the mesh structure.\n"));
430
 
 
431
 
}
432
 
 
433
 
/*----------------------------------------------------------------------------*/
434
 
 
435
 
END_C_DECLS