1 /*============================================================================
2  * \file measures set management.
3  *============================================================================*/
4 
5 /*
6   This file is part of Code_Saturne, a general-purpose CFD tool.
7 
8   Copyright (C) 1998-2021 EDF S.A.
9 
10   This program is free software; you can redistribute it and/or modify it under
11   the terms of the GNU General Public License as published by the Free Software
12   Foundation; either version 2 of the License, or (at your option) any later
13   version.
14 
15   This program is distributed in the hope that it will be useful, but WITHOUT
16   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17   FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
18   details.
19 
20   You should have received a copy of the GNU General Public License along with
21   this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
22   Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 */
24 
25 /*----------------------------------------------------------------------------*/
26 
27 #include "cs_defs.h"
28 
29 /*----------------------------------------------------------------------------
30  * Standard C library headers
31  *----------------------------------------------------------------------------*/
32 
33 #include <assert.h>
34 #include <math.h>
35 #include <stdio.h>
36 #include <stdlib.h>
37 #include <string.h>
38 #include <stdarg.h>
39 #include <float.h>
40 
41 #if defined(HAVE_MPI)
42 #include <mpi.h>
43 #endif
44 
45 /*----------------------------------------------------------------------------
46  * Local headers
47  *----------------------------------------------------------------------------*/
48 
49 #include "bft_mem_usage.h"
50 #include "bft_error.h"
51 #include "bft_printf.h"
52 #include "bft_mem.h"
53 
54 #include "fvm_nodal.h"
55 #include "fvm_point_location.h"
56 
57 #include "cs_base.h"
58 #include "cs_selector.h"
59 #include "cs_mesh_connect.h"
60 #include "cs_post.h"
61 #include "cs_prototypes.h"
62 
63 #include "cs_log.h"
64 #include "cs_map.h"
65 #include "cs_parall.h"
66 #include "cs_mesh_location.h"
67 
68 #include "cs_mesh.h"
69 #include "cs_measures_util.h"
70 
71 #include "cs_field.h"
72 
73 /*----------------------------------------------------------------------------*/
74 
75 BEGIN_C_DECLS
76 
77 /*! \cond DOXYGEN_SHOULD_SKIP_THIS */
78 
79 /*=============================================================================
80  * Local Macro Definitions
81  *============================================================================*/
82 
83 /*============================================================================
84  * Static global variables
85  *============================================================================*/
86 
87 /* Measures set definitions */
88 
89 static int  _n_measures_sets = 0;
90 static int  _n_measures_sets_max = 0;
91 static cs_measures_set_t  *_measures_sets = NULL;
92 static cs_map_name_to_id_t  *_measures_sets_map = NULL;
93 
94 static int  _n_grids = 0;
95 static int  _n_grids_max = 0;
96 static cs_interpol_grid_t  *_grids = NULL;
97 static cs_map_name_to_id_t  *_grids_map = NULL;
98 
99 /*============================================================================
100  * Local Type Definitions
101  *============================================================================*/
102 
103 typedef struct
104 {
105   double val;
106   int    rank;
107 } _cs_base_mpi_double_int_t;
108 
109 /*============================================================================
110  * Private function definitions
111  *============================================================================*/
112 
113 /*----------------------------------------------------------------------------
114  * Create mesh <-> interpolation grid connectivity
115  *
116  * parameters:
117  *   ig             <-- pointer to the interpolation grid structure
118  *----------------------------------------------------------------------------*/
119 
120 static void
_mesh_interpol_create_connect(cs_interpol_grid_t * ig)121 _mesh_interpol_create_connect(cs_interpol_grid_t   *ig)
122 {
123   cs_lnum_t  ii;
124   cs_lnum_t nb_points = ig->nb_points;
125   cs_real_t *coords = ig->coords;
126 
127   cs_lnum_t *location = NULL;
128   float *distance = NULL;
129   fvm_nodal_t *nodal_mesh = NULL;
130   const cs_mesh_t *mesh = cs_glob_mesh;
131 
132 #if defined(HAVE_MPI)
133   _cs_base_mpi_double_int_t val_min, val_in;
134 #endif
135 
136   nodal_mesh = cs_mesh_connect_cells_to_nodal(mesh,
137                                               "temporary",
138                                               false,
139                                               mesh->n_cells,
140                                               NULL);
141 
142   BFT_MALLOC(location, nb_points, cs_lnum_t);
143   BFT_MALLOC(distance, nb_points, float);
144 
145 #   pragma omp parallel for
146   for (ii = 0; ii < nb_points; ii++) {
147     location[ii] = -1;
148     distance[ii] = -1.0;
149   }
150 
151   fvm_point_location_nodal(nodal_mesh,
152                            0,
153                            0.1,
154                            0,
155                            nb_points,
156                            NULL,
157                            coords,
158                            location,
159                            distance);
160 
161 
162 #if defined(HAVE_MPI)
163   if (cs_glob_n_ranks > 1) {
164     for (ii = 0; ii < nb_points; ii++) {
165       if (location[ii] > 0)
166         val_in.val = distance[ii];
167       else
168         val_in.val = DBL_MAX;
169 
170       val_in.rank = cs_glob_rank_id;
171 
172       MPI_Reduce(&val_in, &val_min, 1, MPI_DOUBLE_INT, MPI_MINLOC, 0,
173                   cs_glob_mpi_comm);
174       MPI_Bcast(&val_min.rank, 1, MPI_INT, 0, cs_glob_mpi_comm);
175       MPI_Bcast(&location[ii], 1, MPI_INT, val_min.rank,
176                 cs_glob_mpi_comm);
177 
178       ig->rank_connect[ii] = val_min.rank;
179     }
180   }
181 #endif
182 
183   /* copy location and switch to 0-based */
184 #   pragma omp parallel for
185   for (ii = 0; ii < nb_points; ii++) {
186     ig->cell_connect[ii] = location[ii] - 1;
187   }
188 
189   nodal_mesh = fvm_nodal_destroy(nodal_mesh);
190   BFT_FREE(location);
191   BFT_FREE(distance);
192 }
193 
194 /*----------------------------------------------------------------------------
195  * Interpolate mesh field on interpol grid structure.
196  *
197  * parameters:
198  *   ig                   <-- pointer to the interpolation grid structure
199  *   values_to_interpol   <-- field on mesh (size = n_cells)
200  *   interpolated_values  --> interpolated values on the interpolation grid
201  *                            structure (size = ig->nb_point)
202  *----------------------------------------------------------------------------*/
203 
204 void
cs_interpol_field_on_grid(cs_interpol_grid_t * ig,const cs_real_t * values_to_interpol,cs_real_t * interpoled_values)205 cs_interpol_field_on_grid(cs_interpol_grid_t         *ig,
206                           const cs_real_t            *values_to_interpol,
207                           cs_real_t                  *interpoled_values)
208 {
209   cs_lnum_t ii, jj;
210   int ms_dim = 1;
211   cs_lnum_t nb_points = ig->nb_points;
212   cs_mesh_t *mesh = cs_glob_mesh;
213 
214 #   pragma omp parallel for private(jj)
215   for (ii = 0; ii < nb_points; ii++) {
216     if (ig->cell_connect[ii] > -1 && ig->cell_connect[ii] < mesh->n_cells)
217       for (jj = 0; jj < ms_dim; jj++)
218         interpoled_values[ii*ms_dim +jj] =
219           values_to_interpol[(ig->cell_connect[ii])*ms_dim + jj];
220   }
221 
222 #if defined(HAVE_MPI)
223   if (cs_glob_n_ranks > 1)
224     for (ii = 0; ii < nb_points; ii++)
225       for (jj = 0; jj < ms_dim; jj++)
226         MPI_Bcast(&(interpoled_values[ii*ms_dim + jj]), 1,
227                   CS_MPI_REAL, ig->rank_connect[ii],
228                   cs_glob_mpi_comm);
229 
230 #endif
231 }
232 
233 /*----------------------------------------------------------------------------
234  * Compute a Cressman interpolation on the global mesh.
235  *
236  * parameters:
237  *   ms                   <-- pointer to the measures set structure
238  *                            (values to interpolate)
239  *   interpolated_values  --> interpolated values on the global mesh
240  *                            (size = n_cells or nb_faces)
241  *   id_type              <-- parameter:
242  *                              1: interpolation on volumes
243  *                              2: interpolation on boundary faces
244  *----------------------------------------------------------------------------*/
245 
246 void
cs_cressman_interpol(cs_measures_set_t * ms,cs_real_t * interpolated_values,int id_type)247 cs_cressman_interpol(cs_measures_set_t         *ms,
248                      cs_real_t                 *interpolated_values,
249                      int                        id_type)
250 {
251   cs_lnum_t n_elts = 0;
252   cs_real_t *xyz_cen = NULL;
253   const cs_mesh_t *mesh = cs_glob_mesh;
254   const cs_mesh_quantities_t *mesh_quantities = cs_glob_mesh_quantities;
255 
256   if (id_type == 1) {
257     n_elts = mesh->n_cells;
258     xyz_cen = mesh_quantities->cell_cen;
259   }
260   else if (id_type == 2) {
261     n_elts = mesh->n_b_faces;
262     xyz_cen = mesh_quantities->b_face_cog;
263   }
264   assert(xyz_cen != NULL);
265 
266 # pragma omp parallel for
267   for (cs_lnum_t ii = 0; ii < n_elts; ii++) {
268     cs_real_t total_weight = 0.;
269     cs_real_t interpolated_value = 0.;
270     for (cs_lnum_t jj = 0; jj < ms->nb_measures; jj++) {
271       if (ms->is_cressman[jj] == 1) {
272         cs_real_t dist_x = (xyz_cen[ii*3   ] - ms->coords[jj*3   ])
273                           *ms->inf_radius[jj*3   ];
274         cs_real_t dist_y = (xyz_cen[ii*3 +1] - ms->coords[jj*3 +1])
275                           *ms->inf_radius[jj*3 +1];
276         cs_real_t dist_z = (xyz_cen[ii*3 +2] - ms->coords[jj*3 +2])
277                           *ms->inf_radius[jj*3 +2];
278 
279         cs_real_t r2 = dist_x*dist_x + dist_y*dist_y + dist_z*dist_z;
280 
281         cs_real_t weight = 0.;
282         if (r2/4. <= 700.)
283           weight = exp(-r2/4.);
284 
285         total_weight += weight;
286         interpolated_value += (ms->measures[jj])*weight;
287       }
288     }
289 
290     if (total_weight > 0.)
291       interpolated_values[ii] = interpolated_value/total_weight;
292     else
293       interpolated_values[ii] = 0.;
294   }
295 
296 }
297 
298 /*----------------------------------------------------------------------------
299  * Create an interpolation grid descriptor.
300  *
301  * For measures set with a dimension greater than 1, components are interleaved.
302  *
303  * parameters:
304  *   name        <-- grid name
305  *
306  * returns:
307  *   pointer to new interpolation grid.
308  *
309  *----------------------------------------------------------------------------*/
310 
311 cs_interpol_grid_t *
cs_interpol_grid_create(const char * name)312 cs_interpol_grid_create(const char   *name)
313 {
314   bool reall = true;
315   int grid_id = -1;
316 
317   const char *addr_0 = NULL, *addr_1 = NULL;
318 
319   cs_interpol_grid_t *ig =  NULL;
320 
321   /* Initialize if necessary */
322 
323   if (_grids_map == NULL)
324     _grids_map = cs_map_name_to_id_create();
325   else
326     addr_0 = cs_map_name_to_id_reverse(_grids_map, 0);
327 
328   if (strlen(name) == 0)
329     bft_error(__FILE__, __LINE__, 0,
330               _("Defining a interpolation grid requires a name."));
331 
332   /* Find or insert entry in map */
333 
334   grid_id = cs_map_name_to_id(_grids_map, name);
335 
336   /* Move name pointers of previous grid if necessary
337      (i.e. reallocation of map names array) */
338 
339   addr_1 = cs_map_name_to_id_reverse(_grids_map, 0);
340 
341   if (addr_1 != addr_0) {
342     int i;
343     ptrdiff_t addr_shift = addr_1 - addr_0;
344     for (i = 0; i < grid_id; i++)
345       (_grids + i)->name += addr_shift;
346   }
347 
348   if (grid_id == _n_grids) {
349     _n_grids = grid_id + 1;
350     reall = false;
351   }
352 
353   /* Reallocate grids if necessary */
354 
355   if (_n_grids > _n_grids_max) {
356     if (_n_grids_max == 0)
357       _n_grids_max = 8;
358     else
359       _n_grids_max *= 2;
360     BFT_REALLOC(_grids, _n_grids_max, cs_interpol_grid_t);
361   }
362 
363   /* Assign grid */
364 
365   ig = _grids + grid_id;
366 
367   ig->name = cs_map_name_to_id_reverse(_grids_map, grid_id);
368 
369   ig->id = grid_id;
370 
371   ig->nb_points = 0;
372 
373   if (!reall) {
374     ig->coords = NULL;
375     ig->cell_connect = NULL;
376     ig->rank_connect = NULL;
377   }
378 
379   else {
380     BFT_FREE(ig->coords);
381     if (ig->is_connect) {
382       BFT_FREE(ig->cell_connect);
383 #if defined(HAVE_MPI)
384       if (cs_glob_n_ranks > 1)
385         BFT_FREE(ig->rank_connect);
386 #endif
387     }
388   }
389   ig->is_connect = false;
390 
391   return ig;
392 }
393 
394 /*----------------------------------------------------------------------------
395  * Create interpolation grid.
396  *
397  * parameters:
398  *   ig        <-- pointer to an interpol grid structure
399  *   nb_points <-- number of considered points
400  *   coord     <-- coordonates of considered points
401  *----------------------------------------------------------------------------*/
402 
403 void
cs_interpol_grid_init(cs_interpol_grid_t * ig,const cs_lnum_t nb_points,const cs_real_t * coords)404 cs_interpol_grid_init(cs_interpol_grid_t    *ig,
405                       const cs_lnum_t        nb_points,
406                       const cs_real_t       *coords)
407 {
408   cs_lnum_t ii;
409   BFT_MALLOC(ig->cell_connect, nb_points, cs_lnum_t);
410 #if defined(HAVE_MPI)
411   if (cs_glob_n_ranks > 1)
412     BFT_MALLOC(ig->rank_connect, nb_points, int);
413 #endif
414   BFT_MALLOC(ig->coords, 3*nb_points, cs_real_t);
415 
416 #   pragma omp parallel for
417   for (ii = 0; ii < 3*nb_points; ii++) {
418     ig->coords[ii] = coords[ii];
419   }
420 
421   ig->nb_points = nb_points;
422 
423   _mesh_interpol_create_connect(ig);
424 
425   ig->is_connect = true;
426 }
427 
428 /*----------------------------------------------------------------------------
429  * Create a measures set descriptor.
430  *
431  * For measures set with a dimension greater than 1, components are interleaved.
432  *
433  * parameters:
434  *   name        <-- measures set name
435  *   type_flag   <-- mask of field property and category values (not used yet)
436  *   dim         <-- measure set dimension (number of components)
437  *   interleaved <-- if dim > 1, indicate if field is interleaved
438  *
439  * returns:
440  *   pointer to new measures set.
441  *
442  *----------------------------------------------------------------------------*/
443 
444 cs_measures_set_t *
cs_measures_set_create(const char * name,int type_flag,int dim,bool interleaved)445 cs_measures_set_create(const char   *name,
446                        int           type_flag,
447                        int           dim,
448                        bool          interleaved)
449 {
450   bool reall = true;
451   int measures_set_id = -1;
452 
453   const char *addr_0 = NULL, *addr_1 = NULL;
454 
455   cs_measures_set_t *ms =  NULL;
456 
457   /* Initialize if necessary */
458 
459   if (_measures_sets_map == NULL)
460     _measures_sets_map = cs_map_name_to_id_create();
461   else
462     addr_0 = cs_map_name_to_id_reverse(_measures_sets_map, 0);
463 
464   if (strlen(name) == 0)
465     bft_error(__FILE__, __LINE__, 0,
466               _("Defining a measure set requires a name."));
467 
468   /* Find or insert entry in map */
469 
470   measures_set_id = cs_map_name_to_id(_measures_sets_map, name);
471 
472   /* Move name pointers of previous measure set if necessary
473      (i.e. reallocation of map names array) */
474 
475   addr_1 = cs_map_name_to_id_reverse(_measures_sets_map, 0);
476 
477   if (addr_1 != addr_0) {
478     int i;
479     ptrdiff_t addr_shift = addr_1 - addr_0;
480     for (i = 0; i < measures_set_id; i++)
481       (_measures_sets + i)->name += addr_shift;
482   }
483 
484   if (measures_set_id == _n_measures_sets) {
485     _n_measures_sets = measures_set_id + 1;
486     reall = false;
487   }
488 
489   /* Reallocate measures sets if necessary */
490 
491   if (_n_measures_sets > _n_measures_sets_max) {
492     if (_n_measures_sets_max == 0)
493       _n_measures_sets_max = 8;
494     else
495       _n_measures_sets_max *= 2;
496     BFT_REALLOC(_measures_sets, _n_measures_sets_max, cs_measures_set_t);
497   }
498 
499   /* Assign measure set */
500 
501   ms = _measures_sets + measures_set_id;
502 
503   ms->name = cs_map_name_to_id_reverse(_measures_sets_map, measures_set_id);
504 
505   ms->id = measures_set_id;
506   ms->type = type_flag;
507   ms->dim = dim;
508 
509   if (ms->dim > 1)
510     ms->interleaved = interleaved;
511   else
512     ms->interleaved = true;
513 
514   ms->nb_measures = 0;
515   ms->nb_measures_max = 0;
516 
517   if (!reall) {
518     ms->coords = NULL;
519     ms->measures = NULL;
520     ms->is_cressman = NULL;
521     ms->is_interpol = NULL;
522     ms->inf_radius = NULL;
523     ms->comp_ids = NULL;
524   }
525   else {
526     BFT_FREE(ms->coords);
527     BFT_FREE(ms->measures);
528     BFT_FREE(ms->is_cressman);
529     BFT_FREE(ms->is_interpol);
530     BFT_FREE(ms->inf_radius);
531     BFT_FREE(ms->comp_ids);
532   }
533 
534   return ms;
535 }
536 
537 /*----------------------------------------------------------------------------
538  * (re)Allocate and fill in a measures set structure with an array of measures.
539  *
540  * parameters:
541  *   ms               <-- pointer to the measures set
542  *   nb_measures      <-- number of measures
543  *   is_cressman      <-- for each measure cressman interpolation is:
544  *                          0: not used
545  *                          1: used
546  *   is_interpol      <-- for the interpolation on mesh, each measure is
547  *                          0: not taken into account
548  *                          1: taken into account
549  *   measures_coords  <-- measures spaces coordonates
550  *   measures         <-- measures values (associated to coordinates)
551  *   influence_radius <-- influence radius for interpolation (xyz interleaved)
552  *----------------------------------------------------------------------------*/
553 
554 void
cs_measures_set_map_values(cs_measures_set_t * ms,const cs_lnum_t nb_measures,const int * is_cressman,const int * is_interpol,const cs_real_t * measures_coords,const cs_real_t * measures,const cs_real_t * influence_radius)555 cs_measures_set_map_values(cs_measures_set_t       *ms,
556                            const cs_lnum_t          nb_measures,
557                            const int               *is_cressman,
558                            const int               *is_interpol,
559                            const cs_real_t         *measures_coords,
560                            const cs_real_t         *measures,
561                            const cs_real_t         *influence_radius)
562 {
563   cs_lnum_t  ii;
564   int dim = ms->dim;
565 
566   if (nb_measures != ms->nb_measures) {
567     BFT_REALLOC(ms->measures, nb_measures*dim, cs_real_t);
568     BFT_REALLOC(ms->inf_radius, nb_measures*3, cs_real_t);
569     BFT_REALLOC(ms->coords, nb_measures*3, cs_real_t);
570     BFT_REALLOC(ms->is_cressman, nb_measures, int);
571     BFT_REALLOC(ms->is_interpol, nb_measures, int);
572     ms->nb_measures = nb_measures;
573   }
574 
575   if (dim == 1) {
576 #   pragma omp parallel for
577     for (ii = 0; ii < nb_measures; ii++)
578       ms->measures[ii] = measures[ii];
579   }
580   else {
581     if (ms->interleaved) {
582       cs_lnum_t jj;
583 #   pragma omp parallel for private(jj)
584       for (ii = 0; ii < nb_measures; ii++) {
585         for (jj = 0; jj < dim; jj++)
586           ms->measures[ii*dim + jj] = measures[ii*dim + jj];
587       }
588     }
589     else {
590       cs_lnum_t jj;
591 #   pragma omp parallel for private(jj)
592       for (ii = 0; ii < nb_measures; ii++) {
593         for (jj = 0; jj < dim; jj++)
594           ms->measures[ii*dim + jj] = measures[jj*nb_measures + ii];
595       }
596     }
597   }
598 #   pragma omp parallel for
599   for (ii = 0; ii < nb_measures; ii++) {
600     ms->is_interpol[ii] = is_interpol[ii];
601     ms->is_cressman[ii] = is_cressman[ii];
602   }
603 
604   {
605     cs_lnum_t jj;
606 #   pragma omp parallel for private(jj)
607     for (ii = 0; ii < nb_measures; ii++)
608       for (jj = 0; jj < 3; jj++) {
609         ms->coords[ii*3 + jj] = measures_coords[ii*3 + jj];
610         ms->inf_radius[ii*3 +jj] = influence_radius[ii*3 + jj];
611       }
612   }
613 }
614 
615 /*----------------------------------------------------------------------------
616  * Add new measures to an existing measures set (already declared and
617  * allocated).
618  *
619  * parameters:
620  *   ms               <-- pointer to the existing measures set
621  *   nb_measures      <-- number of new measures
622  *   is_cressman      <-- for each new measure cressman interpolation is:
623  *                          0: not used
624  *                          1: used
625  *   is_interpol      <-- for the interpolation on mesh, each new measure is
626  *                          0: not taken into account
627  *                          1: taken into account
628  *   measures_coords  <-- new measures spaces coordonates
629  *   measures         <-- new measures values (associated to coordonates)
630  *   influence_radius <-- influence radius for interpolation (xyz interleaved)
631  *----------------------------------------------------------------------------*/
632 
633 void
cs_measures_set_add_values(cs_measures_set_t * ms,const cs_lnum_t nb_measures,const int * is_cressman,const int * is_interpol,const cs_real_t * measures_coords,const cs_real_t * measures,const cs_real_t * influence_radius)634 cs_measures_set_add_values(cs_measures_set_t       *ms,
635                            const cs_lnum_t          nb_measures,
636                            const int               *is_cressman,
637                            const int               *is_interpol,
638                            const cs_real_t         *measures_coords,
639                            const cs_real_t         *measures,
640                            const cs_real_t         *influence_radius)
641 {
642   cs_lnum_t  ii;
643   int dim = ms->dim;
644   if (ms->nb_measures + nb_measures > ms->nb_measures_max) {
645     ms->nb_measures_max = 2*(ms->nb_measures + nb_measures);
646 
647     BFT_REALLOC(ms->measures, ms->nb_measures_max*dim, cs_real_t);
648     BFT_REALLOC(ms->coords, ms->nb_measures_max*3, cs_real_t);
649     BFT_REALLOC(ms->is_cressman, ms->nb_measures_max, int);
650     BFT_REALLOC(ms->is_interpol, ms->nb_measures_max, int);
651   }
652 
653   if (dim == 1) {
654 #   pragma omp parallel for
655     for (ii = 0; ii < nb_measures; ii++)
656       ms->measures[ms->nb_measures + ii] = measures[ii];
657   }
658   else {
659     if (ms->interleaved) {
660       cs_lnum_t jj;
661 #   pragma omp parallel for private(jj)
662       for (ii = 0; ii < nb_measures; ii++) {
663         for (jj = 0; jj < dim; jj++)
664           ms->measures[(ii + ms->nb_measures)*dim + jj] = measures[ii*dim + jj];
665       }
666     }
667     else {
668       cs_lnum_t jj;
669 #   pragma omp parallel for private(jj)
670       for (ii = 0; ii < nb_measures; ii++) {
671         for (jj = 0; jj < dim; jj++)
672           ms->measures[ii*dim + jj] = measures[ii*nb_measures + jj];
673       }
674     }
675   }
676 #   pragma omp parallel for
677   for (ii = 0; ii < nb_measures; ii++) {
678     ms->is_interpol[ii + ms->nb_measures] = is_interpol[ii];
679     ms->is_cressman[ii + ms->nb_measures] = is_cressman[ii];
680   }
681 
682   {
683     cs_lnum_t jj;
684 #   pragma omp parallel for private(jj)
685     for (ii = 0; ii < nb_measures; ii++)
686       for (jj = 0; jj < 3; jj++) {
687         ms->coords[(ii + ms->nb_measures)*3 + jj] = measures_coords[ii*3 + jj];
688         ms->inf_radius[(ii + ms->nb_measures)*3 + jj] = influence_radius[ii*3 + jj];
689       }
690   }
691 
692   ms->nb_measures += nb_measures;
693 }
694 
695 /*----------------------------------------------------------------------------
696  * Return a pointer to a measures set based on its id.
697  *
698  * This function requires that a measures set of the given id is defined.
699  *
700  * parameters:
701  *  id <-- measures set id
702  *
703  * return:
704  *   pointer to the measures set structure
705  *
706  *----------------------------------------------------------------------------*/
707 
708 cs_measures_set_t  *
cs_measures_set_by_id(int id)709 cs_measures_set_by_id(int  id)
710 {
711   if (id > -1 && id < _n_measures_sets)
712     return _measures_sets + id;
713   else {
714     bft_error(__FILE__, __LINE__, 0,
715               _("Measure set with id %d is not defined."), id);
716     return NULL;
717   }
718 }
719 
720 /*----------------------------------------------------------------------------
721  * Return a pointer to a grid based on its id.
722  *
723  * This function requires that a grid of the given id is defined.
724  *
725  * parameters:
726  *  id <-- grid id
727  *
728  * return:
729  *   pointer to the grid structure
730  *
731  *----------------------------------------------------------------------------*/
732 
733 cs_interpol_grid_t  *
cs_interpol_grid_by_id(int id)734 cs_interpol_grid_by_id(int  id)
735 {
736   if (id > -1 && id < _n_grids)
737     return _grids + id;
738   else {
739     bft_error(__FILE__, __LINE__, 0,
740               _("Interpol grid with id %d is not defined."), id);
741     return NULL;
742   }
743 }
744 
745 /*----------------------------------------------------------------------------
746  * Return a pointer to a measure set based on its name.
747  *
748  * This function requires that a measure set of the given name is defined.
749  *
750  * parameters:
751  * name <-- measure set name
752  *
753  * return:
754  *   pointer to the measures set structure
755  *----------------------------------------------------------------------------*/
756 
757 cs_measures_set_t  *
cs_measures_set_by_name(const char * name)758 cs_measures_set_by_name(const char  *name)
759 {
760   int id = cs_map_name_to_id_try(_measures_sets_map, name);
761 
762   if (id > -1)
763     return _measures_sets + id;
764   else {
765     bft_error(__FILE__, __LINE__, 0,
766               _("Measure set \"%s\" is not defined."), name);
767     return NULL;
768   }
769 }
770 
771 /*----------------------------------------------------------------------------
772  * Return a pointer to a grid based on its name.
773  *
774  * This function requires that a grid of the given name is defined.
775  *
776  * parameters:
777  * name <-- grid name
778  *
779  * return:
780  *   pointer to the grid structure
781  *----------------------------------------------------------------------------*/
782 
783 cs_interpol_grid_t  *
cs_interpol_grid_by_name(const char * name)784 cs_interpol_grid_by_name(const char  *name)
785 {
786   int id = cs_map_name_to_id_try(_grids_map, name);
787 
788   if (id > -1)
789     return _grids + id;
790   else {
791     bft_error(__FILE__, __LINE__, 0,
792               _("Interpol grid \"%s\" is not defined."), name);
793     return NULL;
794   }
795 }
796 
797 
798 /*----------------------------------------------------------------------------
799  * Destroy all defined measures sets.
800  *----------------------------------------------------------------------------*/
801 
802 void
cs_measures_sets_destroy(void)803 cs_measures_sets_destroy(void)
804 {
805   int i;
806 
807   for (i = 0; i < _n_measures_sets; i++) {
808     cs_measures_set_t  *ms = _measures_sets + i;
809     BFT_FREE(ms->measures);
810     BFT_FREE(ms->coords);
811     BFT_FREE(ms->is_interpol);
812     BFT_FREE(ms->is_cressman);
813     BFT_FREE(ms->inf_radius);
814     BFT_FREE(ms->comp_ids);
815   }
816 
817   BFT_FREE(_measures_sets);
818 
819   cs_map_name_to_id_destroy(&_measures_sets_map);
820 
821   _n_measures_sets = 0;
822   _n_measures_sets_max = 0;
823 }
824 
825 /*----------------------------------------------------------------------------
826  * Destroy all defined grids.
827  *----------------------------------------------------------------------------*/
828 
829 void
cs_interpol_grids_destroy(void)830 cs_interpol_grids_destroy(void)
831 {
832   int i;
833 
834   for (i = 0; i < _n_grids; i++) {
835     cs_interpol_grid_t  *ig = _grids + i;
836     BFT_FREE(ig->coords);
837     BFT_FREE(ig->cell_connect);
838     if (cs_glob_n_ranks > 1)
839       BFT_FREE(ig->rank_connect);
840   }
841 
842   BFT_FREE(_grids);
843 
844   cs_map_name_to_id_destroy(&_grids_map);
845 
846   _n_grids = 0;
847   _n_grids_max = 0;
848 }
849 
850 /*! (DOXYGEN_SHOULD_SKIP_THIS) \endcond */
851 
852 /*============================================================================
853  * Public function definitions for Fortran API
854  *============================================================================*/
855 
856 /*----------------------------------------------------------------------------
857  * Define a measures set.
858  *
859  * Fortran interface; use mestcr;
860  *
861  * subroutine mestcr (name, lname, idim, ilved, imeset)
862  * *****************
863  *
864  * character*       name        : <-- : Measure set name
865  * integer          lname       : <-- : Measure set name length
866  * integer          idim        : <-- : Measures set dimension
867  * integer          ilved       : <-- : 0: not intereaved; 1: interleaved
868  * integer          imesset     : --> : id of defined measures set
869  *----------------------------------------------------------------------------*/
870 
CS_PROCF(mestcr,MESTCR)871 void CS_PROCF(mestcr, MESTCR)
872 (
873  const char   *name,
874  const int    *lname,
875  const int    *idim,
876  const int    *ilved,
877  int          *imeset
878 )
879 {
880   char *bufname;
881   int type_flag = 0;
882   bool interleaved = (*ilved == 0) ? false : true;
883   cs_measures_set_t *ms = NULL;
884 
885   bufname = cs_base_string_f_to_c_create(name, *lname);
886 
887   ms = cs_measures_set_create(bufname,
888                               type_flag,
889                               *idim,
890                               interleaved);
891 
892   cs_base_string_f_to_c_free(&bufname);
893 
894   *imeset = ms->id;
895 }
896 
897 /*----------------------------------------------------------------------------
898  * Define a grid.
899  *
900  * Fortran interface
901  *
902  * subroutine gridcr (name, lname, igrid)
903  * *****************
904  *
905  * character*       name        : <-- : Measure set name
906  * integer          lname       : <-- : Measure set name length
907  * integer          igrid       : --> : id of defined grid
908  *----------------------------------------------------------------------------*/
909 
CS_PROCF(gridcr,GRIDCR)910 void CS_PROCF(gridcr, GRIDCR)
911 (
912  const char     *name,
913  const int      *lname,
914  int            *igrid
915 )
916 {
917   char *bufname;
918   cs_interpol_grid_t *ig = NULL;
919 
920   bufname = cs_base_string_f_to_c_create(name, *lname);
921 
922   ig = cs_interpol_grid_create(bufname);
923 
924   cs_base_string_f_to_c_free(&bufname);
925 
926   *igrid = ig->id;
927 }
928 
929 /*----------------------------------------------------------------------------
930  * (re)Allocate and map values to a measure set.
931  *
932  * Fortran interface
933  *
934  * subroutine mesmap (imeset, inbmes, meset, coords, cressm, interp)
935  * *****************
936  *
937  * integer          imeset      : <-- : Measures set id
938  * integer          inbmes      : <-- : Number of measures
939  * cs_real_t*       meset       : <-- : Pointer to measures values array
940  * cs_real_t*       coords      : <-- : Pointer to measures coordonates array
941  * integer*         cressm      : <-- : Pointer to Cressman interpolation flag
942  * integer*         interp      : <-- : Pointer to interpolation flag
943  * integer*         infrad      : <-- : Influence radius for interpolation
944  *----------------------------------------------------------------------------*/
945 
CS_PROCF(mesmap,MESMAP)946 void CS_PROCF(mesmap, MESMAP)
947 (
948  const int         *imeset,
949  const int         *inbmes,
950  const cs_real_t   *meset,
951  const cs_real_t   *coords,
952  const int         *cressm,
953  const int         *interp,
954  const cs_real_t   *infrad
955 )
956 {
957 
958   cs_measures_set_t *ms = cs_measures_set_by_id(*imeset);
959 
960   cs_measures_set_map_values(ms,
961                              *inbmes,
962                              cressm,
963                              interp,
964                              coords,
965                              meset,
966                              infrad);
967 
968 }
969 
970 /*----------------------------------------------------------------------------
971  * Map a grid grid.
972  *
973  * Fortran interface
974  *
975  * subroutine gridmap (name, lname, igrid)
976  * *****************
977  *
978  * integer          igrid       : <-- : Measures set id
979  * integer          inpts       : <-- : Number of measures
980  * cs_real_t*       coords      : <-- : Pointer to measures coordonates array
981  *----------------------------------------------------------------------------*/
982 
CS_PROCF(grimap,GRIMAP)983 void CS_PROCF(grimap, GRIMAP)
984 (
985  const int         *igrid,
986  const int         *inpts,
987  const cs_real_t   *coords
988 )
989 {
990   cs_interpol_grid_t *ig = cs_interpol_grid_by_id(*igrid);
991 
992   cs_interpol_grid_init(ig,
993                         *inpts,
994                         coords);
995 }
996 
997 /*----------------------------------------------------------------------------
998  * Add values to a measure set.
999  *
1000  * Fortran interface
1001  *
1002  * subroutine mesadd (imeset, inbmes, meset, coords, cressm, interp)
1003  * *****************
1004  *
1005  * integer          imeset      : <-- : Measures set id
1006  * integer          inbmes      : <-- : Number of measures to add
1007  * cs_real_t*       meset       : <-- : Pointer to measures values array
1008  * cs_real_t*       coords      : <-- : Pointer to measures coordonates array
1009  * integer*         cressm      : <-- : Pointer to Cressman interpolation flag
1010  * integer*         interp      : <-- : Pointer to interpolation flag
1011  * integer*         infrad      : <-- : Influence radius for interpolation
1012  *----------------------------------------------------------------------------*/
1013 
CS_PROCF(mesadd,MESADD)1014 void CS_PROCF(mesadd, MESADD)
1015 (
1016  const int         *imeset,
1017  const int         *inbmes,
1018  const cs_real_t   *meset,
1019  const cs_real_t   *coords,
1020  const int         *cressm,
1021  const int         *interp,
1022  const cs_real_t   *infrad
1023 )
1024 {
1025   cs_measures_set_t *ms = cs_measures_set_by_id(*imeset);
1026 
1027   cs_measures_set_add_values(ms,
1028                              *inbmes,
1029                              cressm,
1030                              interp,
1031                              coords,
1032                              meset,
1033                              infrad);
1034 
1035 }
1036 
1037 /*----------------------------------------------------------------------------
1038  * Compute a Cressman interpolation on the global mesh.
1039  *
1040  * Fortran interface
1041  *
1042  * subroutine mscrss (imeset, type, pldval)
1043  * *****************
1044  *
1045  * integer          imeset      : <-- : Measures set id
1046  * integer          type        : <-- : Parameter:
1047  *                                        1: interpolation on volumes
1048  *                                        2: interpolation on boundary faces
1049  * cs_real_t*       pldval      : --> : Interpolated values on the global mesh
1050  *----------------------------------------------------------------------------*/
1051 
CS_PROCF(mscrss,MSCRSS)1052 void CS_PROCF(mscrss, MSCRSS)
1053 (
1054  const int         *imeset,
1055  const int         *type,
1056  cs_real_t         *pldval
1057 )
1058 {
1059   cs_measures_set_t *ms = cs_measures_set_by_id(*imeset);
1060 
1061   cs_cressman_interpol(ms,
1062                        pldval,
1063                        *type);
1064 }
1065 
1066 /*----------------------------------------------------------------------------
1067  * Interpolate calculed field on a grid.
1068  *
1069  * Fortran interface
1070  *
1071  * subroutine gripol (igrid, inval, pldval)
1072  * *****************
1073  *
1074  * integer          igrid       : <-- : Measures set id
1075  * cs_real_t*       inval       : <-- : Values to interpolate
1076  * cs_real_t*       pldval      : --> : Interpolated values on the grid
1077  *----------------------------------------------------------------------------*/
1078 
CS_PROCF(gripol,GRIPOL)1079 void CS_PROCF(gripol, GRIPOL)
1080 (
1081  const int         *igrid,
1082  const cs_real_t   *inval,
1083  cs_real_t         *pldval
1084 )
1085 {
1086   cs_interpol_grid_t *ig = cs_interpol_grid_by_id(*igrid);
1087 
1088   cs_interpol_field_on_grid(ig,
1089                             inval,
1090                             pldval);
1091 }
1092 
1093 /*----------------------------------------------------------------------------
1094  * Destroy measures sets.
1095  *
1096  * Fortran interface
1097  *
1098  * subroutine mestde (void)
1099  * *****************
1100  *----------------------------------------------------------------------------*/
1101 
CS_PROCF(mestde,MESTDE)1102 void CS_PROCF(mestde, MESTDE)(void)
1103 {
1104   cs_measures_sets_destroy();
1105 }
1106 
1107 /*----------------------------------------------------------------------------
1108  * Destroy grids.
1109  *
1110  * Fortran interface
1111  *
1112  * subroutine grides (void)
1113  * *****************
1114  *----------------------------------------------------------------------------*/
1115 
CS_PROCF(grides,GRIDES)1116 void CS_PROCF(grides, GRIDES)(void)
1117 {
1118   cs_interpol_grids_destroy();
1119 }
1120 
1121 /*----------------------------------------------------------------------------*/
1122 
1123 END_C_DECLS
1124