1 /*============================================================================
2 * Parallel interpolation using ParaMEDMEM OverlapDEC
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 <stdarg.h>
34 #include <stdio.h>
35 #include <stdlib.h>
36 #include <string.h>
37 #include <assert.h>
38 #include <math.h>
39
40 #if defined(HAVE_MPI)
41 #include <mpi.h>
42 #endif
43
44 /*----------------------------------------------------------------------------
45 * Local headers
46 *----------------------------------------------------------------------------*/
47
48 #include "bft_error.h"
49 #include "bft_mem.h"
50 #include "bft_printf.h"
51
52 #include "cs_mesh.h"
53 #include "cs_mesh_connect.h"
54 #include "cs_parall.h"
55 #include "cs_prototypes.h"
56 #include "cs_selector.h"
57 #include "cs_timer.h"
58
59 #include "fvm_defs.h"
60 #include "fvm_nodal_from_desc.h"
61
62 /*----------------------------------------------------------------------------
63 * Header for the current file
64 *----------------------------------------------------------------------------*/
65
66 #include "cs_medcoupling_utils.hxx"
67 #include "cs_paramedmem_remapper.h"
68
69 /*----------------------------------------------------------------------------
70 * MEDCOUPLING library headers
71 *----------------------------------------------------------------------------*/
72
73 #if defined(HAVE_PARAMEDMEM) && defined(HAVE_MEDCOUPLING_LOADER)
74
75 #include <MEDLoader.hxx>
76 #include <MEDCoupling_version.h>
77 #include <MEDFileField.hxx>
78 #include <MEDFileFieldMultiTS.hxx>
79
80 #include <MEDFileMesh.hxx>
81 #include <MEDCouplingUMesh.hxx>
82 #include <MEDCouplingField.hxx>
83 #include <MEDCouplingFieldDouble.hxx>
84
85 #include <ParaFIELD.hxx>
86 #include <ParaMESH.hxx>
87 #include <OverlapDEC.hxx>
88 #include <ParaMEDFileMesh.hxx>
89
90 using namespace MEDCoupling;
91
92 #endif
93
94 /*----------------------------------------------------------------------------*/
95
96 /*=============================================================================
97 * Local Structure Definitions
98 *============================================================================*/
99
100 /*----------------------------------------------------------------------------
101 * MEDCoupling/ParaMEDMEM parallel interpolation structure
102 *----------------------------------------------------------------------------*/
103
104 #if defined(HAVE_PARAMEDMEM) && defined(HAVE_MEDCOUPLING_LOADER)
105
106 struct _cs_paramedmem_remapper_t {
107
108 char *name; /* Coupling name */
109
110 cs_medcoupling_mesh_t *local_mesh; /* Local cs_mesh in MED format */
111
112 #if defined(HAVE_PARAMEDMEM) && defined(HAVE_MEDCOUPLING_LOADER)
113 MEDFileUMesh *src_mesh;
114 MEDFileFields *MEDFields;
115 #else
116 void *src_mesh;
117 void *MEDFields;
118 #endif
119
120 int ntsteps;
121 int *iter;
122 int *order;
123 cs_real_t *time_steps;
124
125 /* Bounding sphere for localization */
126 cs_real_t _sphere_cen[3] = {0., 0., 0.};
127 cs_real_t _sphere_rad = 1.e20;
128
129 #if defined(HAVE_PARAMEDMEM) && defined(HAVE_MEDCOUPLING_LOADER)
130 OverlapDEC *odec; /* Overlap data exchange channel */
131 #else
132 void *odec;
133 #endif
134
135 int synced;
136 };
137
138 struct _mesh_transformation_t {
139
140 int type = -1; /* 0: rotation, 1: translation */
141
142 cs_real_t vector[3] = {0., 0., 0.};
143 cs_real_t center[3] = {0., 0., 0.};
144 cs_real_t angle = 0.;
145 };
146
147 static int _n_remappers = 0;
148 static cs_paramedmem_remapper_t **_remapper = NULL;
149
150 static int _n_transformations = 0;
151 static _mesh_transformation_t **_transformations = NULL;
152
153 static bool _transformations_applied = false;
154
155 #endif
156
157 /*============================================================================
158 * Private function definitions
159 *============================================================================*/
160
161 #if defined(HAVE_PARAMEDMEM) && defined(HAVE_MEDCOUPLING_LOADER)
162
163 /*----------------------------------------------------------------------------*/
164 /*!
165 * \brief Create a mesh transformation (rotation or translation)
166 *
167 * \return a mesh_transformation pointer
168 */
169 /*----------------------------------------------------------------------------*/
170
171 static _mesh_transformation_t *
_cs_paramedmem_create_transformation(int type,const cs_real_t center[3],const cs_real_t vector[3],const cs_real_t angle)172 _cs_paramedmem_create_transformation(int type,
173 const cs_real_t center[3],
174 const cs_real_t vector[3],
175 const cs_real_t angle)
176 {
177
178 _mesh_transformation_t *mt = NULL;
179
180 BFT_MALLOC(mt, 1, _mesh_transformation_t);
181
182 mt->type = type;
183 mt->angle = angle;
184
185 for (int i = 0; i < 3; i++) {
186 mt->center[i] = center[i];
187 mt->vector[i] = vector[i];
188 }
189
190 return mt;
191 }
192
193 /*----------------------------------------------------------------------------*/
194 /*!
195 * \brief Reset all mesh transformations
196 *
197 */
198 /*----------------------------------------------------------------------------*/
199
200 static void
_cs_paramedmem_reset_transformations(void)201 _cs_paramedmem_reset_transformations(void)
202 {
203 if (_transformations != NULL) {
204 for (int i = 0; i < _n_transformations; i++) {
205 BFT_FREE(_transformations[i]);
206 }
207 BFT_FREE(_transformations);
208 }
209
210 _n_transformations = 0;
211
212 _transformations_applied = false;
213 }
214
215 /*----------------------------------------------------------------------------*/
216 /*!
217 * \brief Apply the different mesh transformations and update the bounding
218 * sphere
219 */
220 /*----------------------------------------------------------------------------*/
221
222 static void
_cs_paramedmem_apply_transformations(MEDCouplingFieldDouble * field,cs_paramedmem_remapper_t * r)223 _cs_paramedmem_apply_transformations(MEDCouplingFieldDouble *field,
224 cs_paramedmem_remapper_t *r)
225 {
226 if (_transformations_applied == false) {
227 for (int i = 0; i < _n_transformations; i++) {
228 _mesh_transformation_t *mt = _transformations[i];
229 if (mt->type == 0) {
230 /* Rotation */
231 field->getMesh()->rotate(mt->center, mt->vector, mt->angle);
232
233 /* Bounding sphere */
234 cs_real_t c = cos(mt->angle);
235 cs_real_t s = sin(mt->angle);
236
237 cs_real_t norm = cs_math_3_norm(mt->vector);
238 cs_real_t vec[3] = {0., 0., 0.};
239 for (int id = 0; id < 3; id++)
240 vec[id] = mt->vector[id]/norm;
241
242 cs_real_t xyz[3] = {0., 0., 0.};
243 for (int id = 0; id < 3; id++)
244 xyz[id] = r->_sphere_cen[id] - mt->center[id];
245
246 r->_sphere_cen[0] = (vec[0]*vec[0]*(1. - c) + c) * xyz[0]
247 + (vec[0]*vec[1]*(1. - c) - vec[2]*s) * xyz[1]
248 + (vec[0]*vec[2]*(1. - c) + vec[1]*s) * xyz[2]
249 + mt->center[0];
250
251 r->_sphere_cen[1] = (vec[0]*vec[1]*(1. - c) + vec[2]*s) * xyz[0]
252 + (vec[1]*vec[1]*(1. - c) + c) * xyz[1]
253 + (vec[1]*vec[2]*(1. - c) - vec[0]*s) * xyz[2]
254 + mt->center[1];
255
256 r->_sphere_cen[2] = (vec[2]*vec[0]*(1. - c) - vec[1]*s) * xyz[0]
257 + (vec[2]*vec[1]*(1. - c) + vec[0]*s) * xyz[1]
258 + (vec[2]*vec[2]*(1. - c) + c) * xyz[2]
259 + mt->center[2];
260
261 } else if (mt->type == 1) {
262 /* Translation */
263 field->getMesh()->translate(mt->vector);
264
265 /* Bounding sphere */
266 for (int id = 0; id < 3; id++)
267 r->_sphere_cen[id] += mt->vector[id];
268 }
269 }
270 }
271
272 _transformations_applied = true;
273 }
274
275 /*----------------------------------------------------------------------------*/
276 /*!
277 * \brief Returns an array containing the ranks of Code_Saturne processes in
278 * MPI_COMM_WORLD.
279 *
280 * \return array of ranks in MPI_COMM_WORLD
281 */
282 /*----------------------------------------------------------------------------*/
283
284 static int *
_cs_paramedmem_get_mpi_comm_world_ranks(void)285 _cs_paramedmem_get_mpi_comm_world_ranks(void)
286 {
287 /* Global rank of current rank */
288 int my_rank;
289 MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
290
291 /* Size of the local communicator */
292 int mycomm_size;
293 MPI_Comm_size(cs_glob_mpi_comm, &mycomm_size);
294
295 int *world_ranks;
296 BFT_MALLOC(world_ranks, mycomm_size, int);
297
298 MPI_Allgather(&my_rank, 1, MPI_INT, world_ranks, 1, MPI_INT, cs_glob_mpi_comm);
299
300 return world_ranks;
301 }
302
303 /*----------------------------------------------------------------------------*/
304 /*!
305 * \brief Create the parallel remapper structure based on ParaMEDMEM OverlapDEC
306 *
307 * \param[in] name name of the remapper
308 *
309 * \return pointer to cs_paramedmem_remapper_t struct
310 */
311 /*----------------------------------------------------------------------------*/
312
313 static cs_paramedmem_remapper_t *
_cs_paramedmem_overlap_create(const char * name)314 _cs_paramedmem_overlap_create(const char *name)
315 {
316 cs_paramedmem_remapper_t *r = NULL;
317
318 /* Add corresponding coupling to temporary ICoCo couplings array */
319
320 BFT_MALLOC(r, 1, cs_paramedmem_remapper_t);
321
322 BFT_MALLOC(r->name, strlen(name) + 1, char);
323 strcpy(r->name, name);
324
325 r->iter = NULL;
326 r->order = NULL;
327 r->time_steps = NULL;
328 r->ntsteps = -1;
329
330 r->synced = 0;
331
332 /* Local id's */
333 int *cs_ranks = _cs_paramedmem_get_mpi_comm_world_ranks();
334 int cs_comm_size;
335 MPI_Comm_size(cs_glob_mpi_comm, &cs_comm_size);
336
337 std::set<int> grp_ids;
338 for (int i = 0; i < cs_comm_size; i++) {
339 grp_ids.insert(cs_ranks[i]);
340 }
341
342 /* Create the Overlap DEC */
343 r->odec = new OverlapDEC(grp_ids);
344
345 return r;
346 }
347
348 /*----------------------------------------------------------------------------*/
349 /*!
350 * \brief Set the target mesh for interpolation
351 *
352 * \param[in] r pointer to cs_paramedmem_remapper_t struct
353 * \param[in] name mesh name
354 * \param[in] select_criteria selection criteria for needed cells
355 */
356 /*----------------------------------------------------------------------------*/
357
358 static void
_cs_paramedmem_remapper_target_mesh(cs_paramedmem_remapper_t * r,const char * name,const char * select_criteria)359 _cs_paramedmem_remapper_target_mesh(cs_paramedmem_remapper_t *r,
360 const char *name,
361 const char *select_criteria)
362 {
363
364 assert(r != NULL);
365
366 /* Initialization: Only volumes are possible for this option */
367 r->local_mesh = cs_medcoupling_mesh_create(name,
368 select_criteria,
369 3);
370
371 cs_mesh_t *parent_mesh = cs_glob_mesh;
372
373 /* Building the MED representation of the internal mesh */
374 cs_medcoupling_mesh_copy_from_base(parent_mesh, r->local_mesh, 0);
375
376 return;
377
378 }
379
380 /*----------------------------------------------------------------------------*/
381 /*!
382 * \brief Load the mesh parts on each process
383 *
384 * \param[in] r pointer to cs_paramedmem_remapper_t struct
385 * \param[in] file_name name of med file containing the data
386 * \param[in] mesh_name name of the mesh to read in the med file
387 *
388 */
389 /*----------------------------------------------------------------------------*/
390
391 static void
_cs_paramedmem_load_paramesh(cs_paramedmem_remapper_t * r,char * file_name,char * mesh_name)392 _cs_paramedmem_load_paramesh(cs_paramedmem_remapper_t *r,
393 char *file_name,
394 char *mesh_name)
395 {
396 int myPart;
397 MPI_Comm_rank(cs_glob_mpi_comm, &myPart);
398 int nParts;
399 MPI_Comm_size(cs_glob_mpi_comm, &nParts);
400
401 const std::string fname = file_name;
402
403 if (mesh_name != NULL) {
404 const std::string mname = mesh_name;
405
406 // Mesh is stored with -1, -1 indices in MED files
407 r->src_mesh = ParaMEDFileUMesh::New(myPart,
408 nParts,
409 fname,
410 mname,
411 -1,
412 -1);
413 } else {
414 MEDFileMeshes *MeshList = ParaMEDFileMeshes::New(myPart, nParts, fname);
415 const std::string mname = MeshList->getMeshAtPos(0)->getName();
416
417
418 // Mesh is stored with -1, -1 indices in MED files
419 r->src_mesh = ParaMEDFileUMesh::New(myPart,
420 nParts,
421 fname,
422 mname,
423 -1,
424 -1);
425
426 }
427
428 r->src_mesh->forceComputationOfParts();
429
430 MCAuto<MEDFileMeshes> ms = MEDFileMeshes::New();
431 ms = MEDFileMeshes::New();
432 ms->pushMesh(r->src_mesh);
433
434 r->MEDFields = MEDFileFields::LoadPartOf(fname, true, ms);
435 r->MEDFields->loadArrays();
436
437 /* Get Time steps info */
438 MCAuto<MEDFileAnyTypeFieldMultiTS> fts = r->MEDFields->getFieldAtPos(0);
439
440 std::vector< std::pair<int,int> > tio = fts->getIterations();
441
442 r->ntsteps = tio.size();
443
444 BFT_MALLOC(r->time_steps, r->ntsteps, cs_real_t);
445 BFT_MALLOC(r->iter, r->ntsteps, int);
446 BFT_MALLOC(r->order, r->ntsteps, int);
447 for (int i = 0; i < r->ntsteps; i++) {
448 int it = tio[i].first;
449 int ord = tio[i].second;
450
451 r->iter[i] = it;
452 r->order[i] = ord;
453
454 r->time_steps[i] = fts->getTimeStep(it,ord)->getTime(it,ord);
455 }
456 }
457
458 /*----------------------------------------------------------------------------*/
459 /*!
460 * \brief Destroy a remapper
461 *
462 * \param[in] r remapper to destroy
463 */
464 /*----------------------------------------------------------------------------*/
465
466 static void
_cs_paramedmem_remapper_destroy(cs_paramedmem_remapper_t * r)467 _cs_paramedmem_remapper_destroy(cs_paramedmem_remapper_t *r)
468 {
469
470 BFT_FREE(r->name);
471 BFT_FREE(r->src_mesh);
472 BFT_FREE(r->MEDFields);
473 BFT_FREE(r->iter);
474 BFT_FREE(r->order);
475 BFT_FREE(r->time_steps);
476 BFT_FREE(r->odec);
477
478 cs_medcoupling_mesh_destroy(r->local_mesh);
479
480 return;
481 }
482
483 #endif /* !HAVE_PARAMEDMEM && !HAVE_MEDCOUPLING_LOADER */
484
485 /*============================================================================
486 * Public C functions
487 *============================================================================*/
488
489 BEGIN_C_DECLS
490
491 /*----------------------------------------------------------------------------*/
492 /*!
493 * \brief Creates a new cs_paramedmem_remapper_t instance
494 *
495 * \param[in] name name of the remapper
496 * \param[in] sel_criteria cells selection criteria
497 * \param[in] file_name med file name
498 * \param[in] mesh_name name of the mesh in the med file
499 * \param[in] center center of bounding sphere
500 * \param[in] radius radius of bounding sphere
501 *
502 * \return cs_paramedmem_remapper_t struct
503 */
504 /*----------------------------------------------------------------------------*/
505
506 cs_paramedmem_remapper_t *
cs_paramedmem_remapper_create(char * name,const char * sel_criteria,char * file_name,char * mesh_name,cs_real_t center[3],cs_real_t radius)507 cs_paramedmem_remapper_create(char *name,
508 const char *sel_criteria,
509 char *file_name,
510 char *mesh_name,
511 cs_real_t center[3],
512 cs_real_t radius)
513 {
514 cs_paramedmem_remapper_t *r = NULL;
515
516 #if !defined(HAVE_PARAMEDMEM) || !defined(HAVE_MEDCOUPLING_LOADER)
517
518 CS_NO_WARN_IF_UNUSED(name);
519 CS_NO_WARN_IF_UNUSED(sel_criteria);
520 CS_NO_WARN_IF_UNUSED(file_name);
521 CS_NO_WARN_IF_UNUSED(mesh_name);
522 CS_NO_WARN_IF_UNUSED(center);
523 CS_NO_WARN_IF_UNUSED(radius);
524
525 bft_error(__FILE__, __LINE__, 0,
526 _("This function cannot be called without "
527 "MEDCoupling MPI support.\n"));
528
529 #else
530
531 if (_remapper == NULL)
532 BFT_MALLOC(_remapper, 1, cs_paramedmem_remapper_t *);
533 else
534 BFT_REALLOC(_remapper, _n_remappers+1, cs_paramedmem_remapper_t *);
535
536 r = _cs_paramedmem_overlap_create(name);
537
538 _cs_paramedmem_remapper_target_mesh(r, name, sel_criteria);
539
540 _cs_paramedmem_load_paramesh(r, file_name, mesh_name);
541
542 r->_sphere_rad = radius;
543 for (int i = 0; i < 3; i++)
544 r->_sphere_cen[i] = center[i];
545
546 _remapper[_n_remappers] = r;
547 _n_remappers++;
548
549 #endif
550
551 return r;
552 }
553
554 /*----------------------------------------------------------------------------*/
555 /*!
556 * \brief get a remapper by its name
557 *
558 * \param[in] name name of the remapper
559 *
560 * \return pointer to cs_paramedmem_remapper_t struct
561 */
562 /*----------------------------------------------------------------------------*/
563
564 cs_paramedmem_remapper_t *
cs_paramedmem_remapper_by_name_try(const char * name)565 cs_paramedmem_remapper_by_name_try(const char *name)
566 {
567 #if !defined(HAVE_PARAMEDMEM) || !defined(HAVE_MEDCOUPLING_LOADER)
568 CS_NO_WARN_IF_UNUSED(name);
569 bft_error(__FILE__, __LINE__, 0,
570 _("This function cannot be called without "
571 "MEDCoupling MPI support.\n"));
572 #else
573
574 if (_n_remappers > 0) {
575 for (int r_id = 0; r_id < _n_remappers; r_id++) {
576 const char *r_name = _remapper[r_id]->name;
577 if (strcmp(r_name, name) == 0) {
578 return _remapper[r_id];
579
580 }
581 }
582 }
583
584 #endif
585
586 return NULL;
587 }
588
589 /*----------------------------------------------------------------------------*/
590 /*!
591 * \brief Remaps a field from the med file to the local mesh for a given time
592 *
593 * \param[in] r pointer to cs_paramedmem_remapper_t struct
594 * \param[in] field_name name of the field to remap from the file
595 * \param[in] default_val default value for unmapped elements
596 * \param[in] dt time value to use from the file
597 * \param[in] it time iteration to use from the file
598 *
599 * \return cs_real_t pointer containing the new values on target mesh
600 */
601 /*----------------------------------------------------------------------------*/
602
603 cs_real_t *
cs_paramedmem_remap_field_one_time(cs_paramedmem_remapper_t * r,char * field_name,cs_real_t default_val,int dt,int it)604 cs_paramedmem_remap_field_one_time(cs_paramedmem_remapper_t *r,
605 char *field_name,
606 cs_real_t default_val,
607 int dt,
608 int it)
609 {
610
611 cs_real_t *new_vals = NULL;
612
613 #if !defined(HAVE_PARAMEDMEM) || !defined(HAVE_MEDCOUPLING_LOADER)
614
615 CS_NO_WARN_IF_UNUSED(r);
616 CS_NO_WARN_IF_UNUSED(field_name);
617 CS_NO_WARN_IF_UNUSED(default_val);
618 CS_NO_WARN_IF_UNUSED(dt);
619 CS_NO_WARN_IF_UNUSED(it);
620
621 bft_error(__FILE__, __LINE__, 0,
622 _("This function cannot be called without "
623 "MEDCoupling MPI support.\n"));
624
625 #else
626 /* Source Field */
627 const std::string fname(field_name);
628
629 MCAuto<MEDFileAnyTypeField1TS> f =
630 r->MEDFields->getFieldWithName(fname)->getTimeStep(dt,it);
631
632 MCAuto<MEDFileField1TS>
633 sf(MEDCoupling::DynamicCast<MEDFileAnyTypeField1TS,MEDFileField1TS>(f));
634
635 MEDCouplingFieldDouble *tmp_field(sf->field(r->src_mesh));
636
637 MEDCouplingFieldDouble *src_field;
638 if (tmp_field->getTypeOfField() == ON_NODES)
639 src_field = tmp_field->nodeToCellDiscretization();
640 else
641 src_field = tmp_field;
642
643 src_field->setNature(IntensiveMaximum);
644
645 _cs_paramedmem_apply_transformations(src_field, r);
646
647 r->odec->attachSourceLocalField(src_field);
648
649 /* Target Field */
650 MEDCouplingFieldDouble *trg_field
651 = MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME);
652 trg_field->setMesh(r->local_mesh->med_mesh);
653
654 DataArrayDouble *arr = DataArrayDouble::New();
655
656 int fdim = src_field->getNumberOfComponents();
657 arr->alloc(r->local_mesh->n_elts, fdim);
658 trg_field->setArray(arr);
659 trg_field->getArray()->decrRef();
660 trg_field->setNature(IntensiveMaximum);
661 r->odec->attachTargetLocalField(trg_field);
662
663 r->odec->setDefaultValue(default_val);
664 // Sync the DEC if needed
665 if (r->synced != 1) {
666 r->odec->synchronize();
667 r->synced = 1;
668 }
669
670 r->odec->sendData();
671
672 // Write new values
673 const double *val_ptr = trg_field->getArray()->getConstPointer();
674
675 BFT_MALLOC(new_vals, r->local_mesh->n_elts, cs_real_t);
676
677 const cs_lnum_t *new_connec = r->local_mesh->new_to_old;
678
679 const int npts = trg_field->getNumberOfValues();
680
681 const cs_real_3_t * xyzcen = (cs_real_3_t *)cs_glob_mesh_quantities->cell_cen;
682
683 for (int ii = 0; ii < npts; ii++) {
684 int c_id = new_connec[ii];
685 if (cs_math_3_distance(xyzcen[c_id], r->_sphere_cen) < r->_sphere_rad)
686 new_vals[c_id] = val_ptr[ii];
687 else
688 new_vals[c_id] = default_val;
689 }
690
691 #endif
692
693 return new_vals;
694 }
695
696 /*----------------------------------------------------------------------------*/
697 /*!
698 * \brief Interpolate a given field on the local mesh for a given time
699 *
700 * \param[in] r pointer to cs_paramedmem_remapper_t struct
701 * \param[in] field_name name of the field to remap from the file
702 * \param[in] default_val default value for unmapped elements
703 * \param[in] time_choice Choice of the time interpolation.
704 * 0: Value of field interpolated at t=tval from the
705 * med file.
706 * 1: Returns field values for the first time step in
707 * the file. tval is then ignored.
708 * 2: Returns field values for the last time step in
709 * the file. tval is then ignored.
710 * \param[in] tval requested time instant. If time choice is 0 and
711 * tval outside of the file time bounds, return value
712 * will be at the the first time step (if tval < tmin)
713 * or last time step (if tval > tmax)
714 *
715 * \return cs_real_t pointer containing the new values on target mesh
716 */
717 /*----------------------------------------------------------------------------*/
718
719 cs_real_t *
cs_paramedmem_remap_field(cs_paramedmem_remapper_t * r,char * field_name,cs_real_t default_val,int time_choice,double tval)720 cs_paramedmem_remap_field(cs_paramedmem_remapper_t *r,
721 char *field_name,
722 cs_real_t default_val,
723 int time_choice,
724 double tval)
725 {
726 cs_real_t *new_vals = NULL;
727
728 #if !defined(HAVE_PARAMEDMEM) || !defined(HAVE_MEDCOUPLING_LOADER)
729 CS_NO_WARN_IF_UNUSED(r);
730 CS_NO_WARN_IF_UNUSED(field_name);
731 CS_NO_WARN_IF_UNUSED(default_val);
732 CS_NO_WARN_IF_UNUSED(time_choice);
733 CS_NO_WARN_IF_UNUSED(tval);
734
735 bft_error(__FILE__, __LINE__, 0,
736 _("This function cannot be called without "
737 "MEDCoupling MPI support.\n"));
738 #else
739 if ((time_choice == 0 && tval < r->time_steps[0]) ||
740 time_choice == 1 ||
741 r->ntsteps == 1) {
742 /* First instance */
743 int it = r->iter[0];
744 int order = r->order[0];
745 new_vals = cs_paramedmem_remap_field_one_time(r,
746 field_name,
747 default_val,
748 it,
749 order);
750
751 }
752 else if ( (time_choice == 0 && tval > r->time_steps[r->ntsteps-1])
753 || time_choice == 2) {
754 /* Last instance */
755 int it = r->iter[r->ntsteps-1];
756 int order = r->order[r->ntsteps-1];
757
758 new_vals
759 = cs_paramedmem_remap_field_one_time(r, field_name, default_val, it, order);
760
761 } else if (time_choice == 0) {
762 /* A given time within the file time bounds*/
763 int id1 = -1;
764 int id2 = -1;
765 for (int i = 0; i < r->ntsteps-1; i++) {
766 if (tval > r->time_steps[i] && tval < r->time_steps[i+1]) {
767 id1 = i;
768 id2 = i+1;
769 break;
770 }
771 }
772
773 cs_real_t t1 = r->time_steps[id1];
774 cs_real_t t2 = r->time_steps[id2];
775
776 cs_real_t *vals1 = cs_paramedmem_remap_field_one_time(r,
777 field_name,
778 default_val,
779 r->iter[id1],
780 r->order[id1]);
781
782 cs_real_t *vals2 = cs_paramedmem_remap_field_one_time(r,
783 field_name,
784 default_val,
785 r->iter[id2],
786 r->order[id2]);
787
788 BFT_MALLOC(new_vals, r->local_mesh->n_elts, cs_real_t);
789 for (int c_id = 0; c_id < r->local_mesh->n_elts; c_id++) {
790 new_vals[c_id] = vals1[c_id] +
791 (vals2[c_id]-vals1[c_id]) *
792 (tval - t1) / (t2-t1);
793 }
794
795 }
796
797 r->synced = 0;
798 _cs_paramedmem_reset_transformations();
799
800 #endif
801
802 return new_vals;
803 }
804
805 /*----------------------------------------------------------------------------*/
806 /*!
807 * \brief translate the mesh using a given vector
808 *
809 * \param[in] r pointer to the cs_paramedmem_remapper_t struct
810 * \param[in] translation translation vector
811 */
812 /*----------------------------------------------------------------------------*/
813
814 void
cs_paramedmem_remapper_translate(cs_paramedmem_remapper_t * r,cs_real_t translation[3])815 cs_paramedmem_remapper_translate(cs_paramedmem_remapper_t *r,
816 cs_real_t translation[3])
817 {
818 #if !defined(HAVE_PARAMEDMEM) || !defined(HAVE_MEDCOUPLING_LOADER)
819 CS_NO_WARN_IF_UNUSED(r);
820 CS_NO_WARN_IF_UNUSED(translation);
821 bft_error(__FILE__, __LINE__, 0,
822 _("This function cannot be called without "
823 "MEDCoupling MPI support.\n"));
824 #else
825 BFT_REALLOC(_transformations, _n_transformations+1, _mesh_transformation_t *);
826
827 cs_real_t cen[3] = {0.,0.,0.};
828 _transformations[_n_transformations] =
829 _cs_paramedmem_create_transformation(1, cen, translation, 0.);
830
831 _n_transformations++;
832 #endif
833 }
834
835 /*----------------------------------------------------------------------------*/
836 /*!
837 * \brief Rotate the mesh using a center point, axis and angle
838 *
839 * \param[in] r pointer to the cs_paramedmem_remapper_t struct
840 * \param[in] invariant coordinates of the invariant point
841 * \param[in] axis rotation axis vector
842 * \param[in] angle rotation angle in radians
843 */
844 /*----------------------------------------------------------------------------*/
845
846 void
cs_paramedmem_remapper_rotate(cs_paramedmem_remapper_t * r,cs_real_t invariant[3],cs_real_t axis[3],cs_real_t angle)847 cs_paramedmem_remapper_rotate(cs_paramedmem_remapper_t *r,
848 cs_real_t invariant[3],
849 cs_real_t axis[3],
850 cs_real_t angle)
851 {
852 #if !defined(HAVE_PARAMEDMEM) || !defined(HAVE_MEDCOUPLING_LOADER)
853 CS_NO_WARN_IF_UNUSED(r);
854 CS_NO_WARN_IF_UNUSED(invariant);
855 CS_NO_WARN_IF_UNUSED(axis);
856 CS_NO_WARN_IF_UNUSED(angle);
857
858 bft_error(__FILE__, __LINE__, 0,
859 _("This function cannot be called without "
860 "MEDCoupling MPI support.\n"));
861 #else
862 BFT_REALLOC(_transformations, _n_transformations+1, _mesh_transformation_t *);
863
864 _transformations[_n_transformations] =
865 _cs_paramedmem_create_transformation(0, invariant, axis, angle);
866
867 _n_transformations++;
868 #endif
869
870 }
871
872 /*----------------------------------------------------------------------------*/
873 /*!
874 * \brief Destroy all remappers
875 */
876 /*----------------------------------------------------------------------------*/
877
878 void
cs_paramedmem_remapper_destroy_all(void)879 cs_paramedmem_remapper_destroy_all(void)
880 {
881
882 #if defined(HAVE_PARAMEDMEM) && defined(HAVE_MEDCOUPLING_LOADER)
883 for (int r_id = 0; r_id < _n_remappers; r_id++)
884 _cs_paramedmem_remapper_destroy(_remapper[r_id]);
885 #else
886 bft_error(__FILE__, __LINE__, 0,
887 _("This function cannot be called without "
888 "MEDCoupling MPI support.\n"));
889 #endif
890
891 return;
892
893 }
894
895 /*----------------------------------------------------------------------------*/
896
897 END_C_DECLS
898