1 /*============================================================================
2 * Write a nodal representation associated with a mesh and associated
3 * variables to Catalyst objects
4 *============================================================================*/
5
6 /*
7 This file is part of Code_Saturne, a general-purpose CFD tool.
8
9 Copyright (C) 1998-2021 EDF S.A.
10
11 This program is free software; you can redistribute it and/or modify it under
12 the terms of the GNU General Public License as published by the Free Software
13 Foundation; either version 2 of the License, or (at your option) any later
14 version.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
19 details.
20
21 You should have received a copy of the GNU General Public License along with
22 this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
23 Street, Fifth Floor, Boston, MA 02110-1301, USA.
24 */
25
26 /*----------------------------------------------------------------------------*/
27
28 /* On glibc-based systems, define _GNU_SOURCE so as to enable
29 modification of floating-point error exceptions handling;
30 _GNU_SOURCE must be defined before including any headers, to ensure
31 the correct feature macros are defined first. */
32
33 #if defined(__linux__) || defined(__linux) || defined(linux)
34 #define CS_FPE_TRAP
35 #if !defined(_GNU_SOURCE)
36 #define _GNU_SOURCE
37 #endif
38 #endif
39
40 #include "cs_defs.h"
41
42 #if defined(HAVE_CATALYST)
43
44 /*----------------------------------------------------------------------------
45 * Standard C library headers
46 *----------------------------------------------------------------------------*/
47
48 #include <assert.h>
49 #include <stdio.h>
50 #include <stdlib.h>
51 #include <string.h>
52
53 /*----------------------------------------------------------------------------
54 * Catalyst and VTK library headers
55 *----------------------------------------------------------------------------*/
56
57 #include <vtkCellType.h>
58
59 #include <vtkCellData.h>
60 #include <vtkCharArray.h>
61 #include <vtkFloatArray.h>
62 #include <vtkMultiBlockDataSet.h>
63 #include <vtkIdList.h>
64 #include <vtkIdTypeArray.h>
65 #include <vtkImageData.h>
66 #include <vtkInformation.h>
67 #include <vtkMPI.h>
68 #include <vtkObjectFactory.h>
69 #include <vtkPointData.h>
70 #include <vtkPolyData.h>
71 #include <vtkUnstructuredGrid.h>
72 #include <vtkPVConfig.h> // vtkPVVersion.h preferred from version 5.10 on.
73
74 #include <vtkCPDataDescription.h>
75 #include <vtkCPInputDataDescription.h>
76 #include <vtkCPProcessor.h>
77 #include <vtkCPPythonScriptPipeline.h>
78 #if defined(HAVE_VTKCPPYTHONSCRIPTV2PIPELINE_H)
79 #include <vtkCPPythonScriptV2Pipeline.h>
80 #endif
81 #include <vtkSmartPointer.h>
82
83 #include <vtkXMLUnstructuredGridWriter.h>
84 #include <vtkUnstructuredGridWriter.h>
85 #include <vtkGenericDataObjectWriter.h>
86 #include <vtkDoubleArray.h>
87
88 #include <vtkFileOutputWindow.h>
89
90 /*----------------------------------------------------------------------------
91 * Local headers
92 *----------------------------------------------------------------------------*/
93
94 #include "bft_error.h"
95 #include "bft_mem.h"
96 #include "bft_printf.h"
97
98 #include "fvm_defs.h"
99 #include "fvm_convert_array.h"
100 #include "fvm_io_num.h"
101 #include "fvm_nodal.h"
102 #include "fvm_nodal_priv.h"
103 #include "fvm_writer_priv.h"
104
105 #include "cs_block_dist.h"
106 #include "cs_file.h"
107 #include "cs_parall.h"
108 #include "cs_part_to_block.h"
109
110 /*----------------------------------------------------------------------------
111 * Header for the current file
112 *----------------------------------------------------------------------------*/
113
114 #include "fvm_to_catalyst.h"
115
116 /*=============================================================================
117 * Local Macro Definitions
118 *============================================================================*/
119
120 /* Macro for version test (tests are done on ParaView and not VTK version,
121 as ParaView usually includes its own VTK, and ParaView 5.0 seems to
122 indicate VTK 7.1 just like 5.1, but that version did not contain
123 SetTypedTuple or have SetTupleValue depecated for vtkTypedDataArray) */
124
125 #define CS_PV_VERSION (PARAVIEW_VERSION_MAJOR*10 + PARAVIEW_VERSION_MINOR)
126
127 /*----------------------------------------------------------------------------
128 * Catalyst field structure
129 *----------------------------------------------------------------------------*/
130
131 typedef struct {
132
133 int mesh_id; /* Associated mesh structure id */
134 int dim; /* Field dimension */
135 vtkUnstructuredGrid *f; /* Pointer to VTK writer fields */
136
137 } fvm_catalyst_field_t;
138
139 /*----------------------------------------------------------------------------
140 * Catalyst writer/reader structure
141 *----------------------------------------------------------------------------*/
142
143 typedef struct {
144
145 char *name; /* Writer name */
146
147 int rank; /* Rank of current process
148 in communicator */
149 int n_ranks; /* Size of communicator */
150
151 vtkMultiBlockDataSet *mb; /* Associated dataset */
152
153 fvm_writer_time_dep_t time_dependency; /* Mesh time dependency */
154
155 int n_fields; /* Number of fields */
156 fvm_catalyst_field_t **fields; /* Array of field helper
157 structures */
158
159 int time_step; /* Latest time step */
160 double time_value; /* Latest time value */
161
162 char *input_name; /* input name, or NULL for
163 default */
164 bool private_comm; /* Use private communicator */
165 bool ensight_names; /* Use EnSight rules for
166 field names */
167
168 bool modified; /* Has output been added since
169 last coprocessing ? */
170
171 } fvm_to_catalyst_t;
172
173 /*============================================================================
174 * Static global variables
175 *============================================================================*/
176
177 int _n_writers = 0;
178
179 int _n_scripts = 0; /* Number of scripts */
180 char **_scripts = NULL; /* List of scripts */
181
182 vtkCPProcessor *_processor = NULL; /* Co processor */
183
184 #if defined(HAVE_MPI)
185 MPI_Comm _comm = MPI_COMM_NULL; /* Associated communicator */
186 MPI_Comm _reference_comm = MPI_COMM_NULL; /* Reference communicator */
187 #endif
188
189 /*============================================================================
190 * Private function definitions
191 *============================================================================*/
192
193 /*----------------------------------------------------------------------------
194 * Check if a script is a Catalyst script
195 *
196 * The script only does cursory checks, so may return false positives
197 * in some cases, but using "incorrect" scripts for Catalyst only
198 * leads to extra warnings, which should be ok pending a more general
199 * file checker.
200 *
201 * parameters:
202 * path <-- scripts path
203 *
204 * returns:
205 * 1 if script is a Catalyst V1, 2 if it seems to be a V2 script,
206 * script, 0 otherwise
207 *----------------------------------------------------------------------------*/
208
209 static int
_check_script_is_catalyst(const char * path)210 _check_script_is_catalyst(const char *path)
211 {
212 assert(path != NULL);
213 int retval = 0;
214
215 FILE *fp = fopen(path, "r");
216
217 if (fp == NULL)
218 return retval;
219
220 int checks[] = {false, false, false};
221 int n_checks = 0;
222 bool import_catalyst = false;
223
224 const char *check_strings[]
225 = {"CreateCoProcessor(",
226 "RequestDataDescription(",
227 "DoCoProcessing("};
228
229 /* Note: we could simply check for "from paraview import coprocessing"
230 in most cases but in cases of non-autogenerated scripts, this might
231 appear in a different form, so we limit ourselves to the main
232 required methods */
233
234 while (true) {
235 char buffer[1024];
236 char *e = buffer+1024;
237
238 char *s = fgets(buffer, 1024, fp);
239 if (s == NULL) break;
240
241 while (s < e && *s == ' ' && *s == '\t') /* skip initial whitespace */
242 s++;
243
244 if (strncmp(s, "def", 3) == 0) {
245 s += 3;
246 if (*s == ' ' || *s == '\t') {
247 while (*s == ' ' && *s == '\t' && *s != '\0')
248 s++;
249 /* Remove whitespace */
250 size_t l = strlen(s);
251 size_t i, j;
252 for (i = 0, j = 0 ; i < l ; i++) {
253 if (s[i] != ' ' && s[i] != '\t')
254 s[j++] = s[i];
255 }
256
257 for (i = 0; i < 3; i++) {
258 if (strncmp(s, check_strings[i], strlen(check_strings[i])) == 0) {
259 if (checks[i] == false) {
260 checks[i] = true;
261 n_checks += 1;
262 }
263 }
264 }
265
266 }
267 }
268 else if (strncmp(s, "from", 4) == 0) {
269
270 /* cleanup whitespace */
271 int n_space = 0;
272 int i = 0;
273 for (int j = 0; j < 1024 && s+j < e && s[j] != '\0'; j++) {
274 if (s[j] == ' ' || s[j] == '\t') {
275 if (n_space < 1)
276 s[i++] = ' ';
277 n_space += 1;
278 }
279 else {
280 s[i++] = s[j];
281 n_space = 0;
282 }
283 }
284
285 if (strncmp(s, "from paraview import catalyst", 29) == 0) {
286 import_catalyst = true;
287 }
288
289 }
290
291 if (n_checks == 3) {
292 retval = 1;
293 break;
294 }
295 else if (import_catalyst) {
296 retval = 2;
297 break;
298 }
299
300 /* Check for end of line; if not present, continue reading from buffer */
301
302 while (s < e && *s != '\0' && *s != '\n')
303 s++;
304 while (s >= e) {
305 s = fgets(buffer, 1024, fp);
306 if (s == NULL) break;
307 while (s < e && *s != '\0' && *s != '\n')
308 s++;
309 }
310 }
311
312 fclose(fp);
313
314 return retval;
315 }
316
317 /*----------------------------------------------------------------------------
318 * Add a Catalyst V2 pipeline if not already present.
319 *
320 * parameters:
321 * path <-- V2 pipeline file path
322 *
323 * returns:
324 * id of pipeline script in list, or -1 if not valid
325 *----------------------------------------------------------------------------*/
326
327 static int
_add_v2_pipeline(const char * path)328 _add_v2_pipeline(const char *path)
329 {
330 assert(path != NULL);
331
332 /* Check that we did not already add this file */
333
334 for (int i = 0; i < _n_scripts; i++) {
335 if (strcmp(_scripts[i], path) == 0)
336 return i;
337 }
338
339 #if defined(HAVE_VTKCPPYTHONSCRIPTV2PIPELINE_H)
340
341 /* Create Catalyst pipeline and check the file is valid */
342
343 vtkNew<vtkCPPythonScriptV2Pipeline> pipeline;
344 if (pipeline->Initialize(path) == false) {
345 bft_printf(_("\nFile \"%s\"\n"
346 "does not seem to contain a Catalyst Python pipeline.\n"),
347 path);
348 return -1;
349 }
350
351 int id = _n_scripts;
352 _n_scripts += 1;
353
354 BFT_REALLOC(_scripts, _n_scripts, char *);
355
356 size_t l = strlen(path);
357 BFT_MALLOC(_scripts[id], l+1, char);
358 strncpy(_scripts[id], path, l);
359 _scripts[id][l] = '\0';
360
361 /* pipeline->SetGhostLevel(1); */
362 _processor->AddPipeline(pipeline);
363
364 return id;
365
366 #else
367
368 bft_printf
369 (_("\nFile \"%s\"\n"
370 "seems to contain a Catalyst V2 Python pipeline but the version.\n"
371 "of Catalyst linked with only supports version 1 pipelines.\n\n"
372 "You may need to use a more recent Catalyst version or downgrade\n"
373 "to earlier V1 pipelines.\n"),
374 path);
375
376 return -1;
377
378 #endif /* defined(HAVE_VTKCPPYTHONSCRIPTV2PIPELINE_H) */
379
380 }
381
382 /*----------------------------------------------------------------------------
383 * Add a Catalyst script if not already present.
384 *
385 * parameters:
386 * path <-- scripts path
387 *
388 * returns:
389 * id of script in list, or -1 if not valid
390 *----------------------------------------------------------------------------*/
391
392 static int
_add_script(const char * path)393 _add_script(const char *path)
394 {
395 assert(path != NULL);
396
397 int is_catalyst = 0;
398 int rank = 0, n_ranks = 1;
399
400 #if defined(HAVE_MPI)
401 if (_comm != MPI_COMM_NULL) {
402 MPI_Comm_rank(_comm, &rank);
403 MPI_Comm_size(_comm, &n_ranks);
404 }
405 #endif
406
407 if (rank < 1)
408 is_catalyst = _check_script_is_catalyst(path);
409
410 #if defined(HAVE_MPI)
411 if (n_ranks > 1)
412 MPI_Bcast(&is_catalyst, 1, MPI_INT, 0, _comm);
413 #endif
414
415 if (is_catalyst < 1)
416 return -1;
417
418 if (is_catalyst == 2) {
419 _add_v2_pipeline(path);
420 return 2;
421 }
422
423 for (int i = 0; i < _n_scripts; i++) {
424 if (strcmp(_scripts[i], path) == 0)
425 return i;
426 }
427
428 int id = _n_scripts;
429 _n_scripts += 1;
430
431 BFT_REALLOC(_scripts, _n_scripts, char *);
432
433 size_t l = strlen(path);
434 BFT_MALLOC(_scripts[id], l+1, char);
435 strncpy(_scripts[id], path, l);
436 _scripts[id][l] = '\0';
437
438 /* Add Catalyst pipeline */
439
440 vtkNew<vtkCPPythonScriptPipeline> pipeline;
441 if (pipeline->Initialize(path) == 0)
442 bft_error(__FILE__, __LINE__, 0,
443 _("Error initializing pipeline from \"%s\""), path);
444
445 /* pipeline->SetGhostLevel(1); */
446 _processor->AddPipeline(pipeline);
447
448 return id;
449 }
450
451 /*----------------------------------------------------------------------------
452 * Add Catalyst scripts from directoty if not already present.
453 *
454 * Currently assumes all Python files in the given directory
455 * are Catalyst scripts.
456 *
457 * parameters:
458 * dir_path <-- directory path
459 *----------------------------------------------------------------------------*/
460
461 static void
_add_dir_scripts(const char * dir_path)462 _add_dir_scripts(const char *dir_path)
463 {
464 char **dir_files = cs_file_listdir(dir_path);
465
466 for (int i = 0; dir_files[i] != NULL; i++) {
467
468 const char *file_name = dir_files[i];
469 const char *ext = NULL;
470 int l_ext = 0;
471
472 /* Find extension */
473 for (int j = strlen(file_name) - 1; j > -1; j--) {
474 l_ext++;
475 if (file_name[j] == '.') {
476 ext = file_name + j;
477 break;
478 }
479 }
480 if (ext == NULL)
481 continue;
482
483 /* Filter: Python files only */
484 if (l_ext == 3 && strncmp(ext, ".py", 3) == 0) {
485 char *tmp_name = NULL;
486 BFT_MALLOC(tmp_name,
487 strlen(dir_path) + 1 + strlen(file_name) + 1,
488 char);
489 sprintf(tmp_name, "%s/%s", dir_path, file_name);
490 _add_script(tmp_name);
491 BFT_FREE(tmp_name);
492 }
493
494 #if defined(HAVE_VTKCPPYTHONSCRIPTV2PIPELINE_H)
495
496 /* Filter: Catalyst V2 pipeline might be in ".zip" form */
497 else if (l_ext == 4 && strncmp(ext, ".zip", 4) == 0) {
498 char *tmp_name = NULL;
499 BFT_MALLOC(tmp_name,
500 strlen(dir_path) + 1 + strlen(file_name) + 1,
501 char);
502 sprintf(tmp_name, "%s/%s", dir_path, file_name);
503 _add_v2_pipeline(tmp_name);
504 BFT_FREE(tmp_name);
505 }
506 #endif /* defined(HAVE_VTKCPPYTHONSCRIPTV2PIPELINE_H) */
507
508 BFT_FREE(dir_files[i]);
509 }
510
511 BFT_FREE(dir_files);
512 }
513
514 /*----------------------------------------------------------------------------
515 * Initialize coprocessor.
516 *
517 * parameters:
518 * private_comm <-- if true, use dedicated communicator
519 * comm <-- associated MPI communicator.
520 *----------------------------------------------------------------------------*/
521
522 #if defined(HAVE_MPI)
523 static void
_init_coprocessor(bool private_comm,MPI_Comm comm)524 _init_coprocessor(bool private_comm,
525 MPI_Comm comm)
526 #else
527 _init_coprocessor(void)
528 #endif
529 {
530 int mpi_flag = 0;
531 int mpi_rank = -1;
532
533 #if defined(HAVE_MPI)
534
535 MPI_Initialized(&mpi_flag);
536
537 if (mpi_flag) {
538 if (comm != _reference_comm) {
539 if (comm != MPI_COMM_NULL && _reference_comm != MPI_COMM_NULL)
540 bft_error(__FILE__, __LINE__, 0,
541 _("All Catalyst writers must use the same MPI communicator"));
542 }
543 }
544
545 #endif
546
547 if (_processor == NULL) {
548
549 _processor = vtkCPProcessor::New();
550
551 #if defined(HAVE_MPI)
552
553 if (mpi_flag) {
554
555 _reference_comm = comm;
556 if (private_comm && comm != MPI_COMM_NULL)
557 MPI_Comm_dup(comm, &(_comm));
558 else
559 _comm = comm;
560
561 if (comm != MPI_COMM_NULL)
562 MPI_Comm_rank(_comm, &mpi_rank);
563
564 vtkMPICommunicatorOpaqueComm vtk_comm
565 = vtkMPICommunicatorOpaqueComm(&_comm);
566 _processor->Initialize(vtk_comm);
567
568 }
569
570 #endif
571
572 if (!mpi_flag)
573 _processor->Initialize();
574
575 vtkFileOutputWindow *log_output = vtkFileOutputWindow::New();
576 if (mpi_rank < 1)
577 log_output->SetFileName("./catalyst.log");
578 else
579 log_output->SetFileName("/dev/null");
580 vtkFileOutputWindow::SetInstance(log_output);
581
582 }
583 }
584
585 /*----------------------------------------------------------------------------
586 * Finalize coprocessor.
587 *
588 * parameters:
589 * private_comm <-- if true, use dedicated communicator
590 * comm <-- associated MPI communicator.
591 *----------------------------------------------------------------------------*/
592
593 static void
_free_coprocessor(void)594 _free_coprocessor(void)
595 {
596 if (_processor != NULL && _n_writers < 2) {
597
598 _processor->Finalize();
599
600 /* Workaround for segmentation fault observed on a Red Hat 8.3 system
601 with gcc 8.3.1, at coprocessor destruction, not observed on other
602 machines. */
603
604 bool cp_delete = true;
605
606 if (_n_scripts > 0) {
607 const char *s = getenv("CS_PV_CP_DELETE_CRASH_WORKAROUND");
608 if (s != NULL) {
609 if (atoi(s) > 0)
610 cp_delete = false;
611 }
612 }
613 if (cp_delete)
614 _processor->Delete();
615
616 _processor = NULL;
617
618 for (int i = 0; i < _n_scripts; i++)
619 BFT_FREE(_scripts[i]);
620 BFT_FREE(_scripts);
621
622 #if defined(HAVE_MPI)
623 {
624 if (_comm != _reference_comm && _comm != MPI_COMM_NULL)
625 MPI_Comm_free(&_comm);
626 }
627 #endif
628
629 }
630 }
631
632 /*----------------------------------------------------------------------------
633 * Return the Catalyst mesh id associated with a given mesh name,
634 * or -1 if no association found.
635 *
636 * parameters:
637 * writer <-- writer structure
638 * mesh_name <-- mesh name
639 *
640 * returns:
641 * Catalyst mesh id, or -1 if Catalyst mesh name is not associated
642 * with this writer structure in FVM
643 *----------------------------------------------------------------------------*/
644
645 static int
_get_catalyst_mesh_id(fvm_to_catalyst_t * writer,const char * mesh_name)646 _get_catalyst_mesh_id(fvm_to_catalyst_t *writer,
647 const char *mesh_name)
648 {
649 int i;
650 int retval = -1;
651
652 assert(writer != NULL);
653
654 int nb = writer->mb->GetNumberOfBlocks();
655 for (i = 0; i < nb; i++) {
656 if (strcmp
657 (mesh_name,
658 writer->mb->GetMetaData(i)->Get(vtkCompositeDataSet::NAME())) == 0)
659 break;
660 }
661
662 if (i < nb)
663 retval = i;
664
665 return retval;
666 }
667
668 /*----------------------------------------------------------------------------
669 * Create a Catalyst mesh structure.
670 *
671 * parameters:
672 * writer <-- Catalyst structure.
673 * mesh <-- FVM mesh structure.
674 *
675 * returns:
676 * Catalyst mesh id.
677 *----------------------------------------------------------------------------*/
678
679 static int
_add_catalyst_mesh(fvm_to_catalyst_t * writer,const fvm_nodal_t * mesh)680 _add_catalyst_mesh(fvm_to_catalyst_t *writer,
681 const fvm_nodal_t *mesh)
682 {
683 assert(writer != NULL);
684
685 /* Add a new Catalyst mesh structure */
686
687 int id = writer->mb->GetNumberOfBlocks();
688 vtkSmartPointer<vtkUnstructuredGrid> ugrid
689 = vtkSmartPointer<vtkUnstructuredGrid>::New();
690
691 writer->mb->SetBlock(id, ugrid);
692 writer->mb->GetMetaData(id)->Set(vtkCompositeDataSet::NAME(), mesh->name);
693
694 return id;
695 }
696
697 /*----------------------------------------------------------------------------
698 * Define VTK geometrical element type according to FVM element type
699 *
700 * parameters:
701 * fvm_elt_type <-- pointer to fvm element type.
702 *
703 * return:
704 * med geometrical element type.
705 *----------------------------------------------------------------------------*/
706
707 static VTKCellType
_get_norm_elt_type(const fvm_element_t fvm_elt_type)708 _get_norm_elt_type(const fvm_element_t fvm_elt_type)
709 {
710 VTKCellType norm_elt_type;
711
712 switch (fvm_elt_type) {
713
714 case FVM_EDGE:
715 norm_elt_type = VTK_LINE;
716 break;
717
718 case FVM_FACE_TRIA:
719 norm_elt_type = VTK_TRIANGLE;
720 break;
721
722 case FVM_FACE_QUAD:
723 norm_elt_type = VTK_QUAD;
724 break;
725
726 case FVM_FACE_POLY:
727 norm_elt_type = VTK_POLYGON;
728 break;
729
730 case FVM_CELL_TETRA:
731 norm_elt_type = VTK_TETRA;
732 break;
733
734 case FVM_CELL_PYRAM:
735 norm_elt_type = VTK_PYRAMID;
736 break;
737
738 case FVM_CELL_PRISM:
739 norm_elt_type = VTK_WEDGE;
740 break;
741
742 case FVM_CELL_HEXA:
743 norm_elt_type = VTK_HEXAHEDRON;
744 break;
745
746 case FVM_CELL_POLY:
747 norm_elt_type = VTK_POLYHEDRON;
748 break;
749
750 default:
751 norm_elt_type = VTK_EMPTY_CELL;
752 bft_error(__FILE__, __LINE__, 0,
753 "_get_norm_elt_type(): "
754 "No association with VTK element type has been found\n"
755 "FVM element type: \"%i\"\n",
756 (int)fvm_elt_type);
757
758 } /* End of switch on element type */
759
760 return norm_elt_type;
761 }
762
763 /*----------------------------------------------------------------------------
764 * Get vertex order to describe Catalyst element type.
765 *
766 * parameters:
767 * norm_elt_type <-- Catalyst element type.
768 * vertex_order --> Pointer to vertex order array (0 to n-1).
769 *----------------------------------------------------------------------------*/
770
771 static void
_get_vertex_order(VTKCellType norm_elt_type,int * vertex_order)772 _get_vertex_order(VTKCellType norm_elt_type,
773 int *vertex_order)
774 {
775 switch(norm_elt_type) {
776
777 case VTK_LINE:
778 vertex_order[0] = 0;
779 vertex_order[1] = 1;
780 break;
781
782 case VTK_TRIANGLE:
783 vertex_order[0] = 0;
784 vertex_order[1] = 1;
785 vertex_order[2] = 2;
786 break;
787
788 case VTK_QUAD:
789 vertex_order[0] = 0;
790 vertex_order[1] = 1;
791 vertex_order[2] = 2;
792 vertex_order[3] = 3;
793 break;
794
795 case VTK_TETRA:
796 vertex_order[0] = 0;
797 vertex_order[1] = 1;
798 vertex_order[2] = 2;
799 vertex_order[3] = 3;
800 break;
801
802 case VTK_PYRAMID:
803 vertex_order[0] = 0;
804 vertex_order[1] = 1;
805 vertex_order[2] = 2;
806 vertex_order[3] = 3;
807 vertex_order[4] = 4;
808 break;
809
810 case VTK_WEDGE:
811 vertex_order[0] = 0;
812 vertex_order[1] = 2;
813 vertex_order[2] = 1;
814 vertex_order[3] = 3;
815 vertex_order[4] = 5;
816 vertex_order[5] = 4;
817 break;
818
819 case VTK_HEXAHEDRON:
820 vertex_order[0] = 0;
821 vertex_order[1] = 1;
822 vertex_order[2] = 2;
823 vertex_order[3] = 3;
824 vertex_order[4] = 4;
825 vertex_order[5] = 5;
826 vertex_order[6] = 6;
827 vertex_order[7] = 7;
828 break;
829
830 case VTK_POLYGON:
831 vertex_order[0] = -1;
832 break;
833
834 case VTK_POLYHEDRON:
835 vertex_order[0] = -1;
836 break;
837
838 default:
839 bft_error(__FILE__, __LINE__, 0,
840 "_get_vertex_order(): No associated FVM element type known\n"
841 "VTK element type: \"%i\"\n",
842 (int)norm_elt_type);
843 }
844
845 return;
846 }
847
848 /*----------------------------------------------------------------------------
849 * Return the Catalyst field id associated with given mesh and field names,
850 * or -1 if no association found.
851 *
852 * parameters:
853 * writer <-- Catalyst writer structure.
854 * fieldname <-- input fieldname.
855 * mesh_id <-- id of associated mesh in structure.
856 * dim <-- number of field components.
857 * location <-- mesh location (cells or vertices).
858 * type <-- field type (cells, nodes)
859 * td <-- time discretization type
860 *
861 * returns
862 * field id in writer structure, or -1 if not found
863 *----------------------------------------------------------------------------*/
864
865 static int
_get_catalyst_field_id(fvm_to_catalyst_t * writer,const char * fieldname,int mesh_id,int dim,cs_datatype_t datatype,fvm_writer_var_loc_t location)866 _get_catalyst_field_id(fvm_to_catalyst_t *writer,
867 const char *fieldname,
868 int mesh_id,
869 int dim,
870 cs_datatype_t datatype,
871 fvm_writer_var_loc_t location)
872 {
873 CS_UNUSED(dim);
874 CS_UNUSED(datatype);
875
876 int i;
877
878 for (i = 0; i < writer->n_fields; ++i){
879
880 vtkFieldData *fData_ptr ;
881
882 if (writer->fields[i]->mesh_id != mesh_id)
883 continue;
884
885 if (location == FVM_WRITER_PER_NODE)
886 fData_ptr = (writer->fields[i])->f->GetPointData();
887 else
888 fData_ptr = writer->fields[i]->f->GetCellData();
889
890 if (fData_ptr->HasArray(fieldname)){
891 break;
892 }
893 }
894
895 if (i < writer->n_fields)
896 return i;
897
898 return -1;
899 }
900
901 /*----------------------------------------------------------------------------
902 * Create a Catalyst field structure.
903 *
904 * parameters:
905 * writer <-- Catalyst writer structure.
906 * fieldname <-- input fieldname.
907 * mesh_id <-- id of associated mesh in structure.
908 * dim <-- number of field components.
909 * location <-- mesh location (cells or vertices).
910 *
911 * returns
912 * field id in writer structure
913 *----------------------------------------------------------------------------*/
914
915 static int
_add_catalyst_field(fvm_to_catalyst_t * writer,const char * fieldname,int mesh_id,int dim,cs_datatype_t datatype,fvm_writer_var_loc_t location)916 _add_catalyst_field(fvm_to_catalyst_t *writer,
917 const char *fieldname,
918 int mesh_id,
919 int dim,
920 cs_datatype_t datatype,
921 fvm_writer_var_loc_t location)
922 {
923 CS_UNUSED(datatype);
924
925 int f_id = writer->n_fields;
926
927 BFT_REALLOC(writer->fields,
928 writer->n_fields+ 1,
929 fvm_catalyst_field_t *);
930
931 vtkUnstructuredGrid *f = NULL;
932
933 const int dest_dim = (dim == 6) ? 9 : dim;
934
935 if (writer->mb->GetMetaData(mesh_id) != NULL) {
936
937 f = vtkUnstructuredGrid::SafeDownCast(vtkDataSet::SafeDownCast
938 (writer->mb->GetBlock(mesh_id)));
939
940 vtkNew<vtkDoubleArray> tmp;
941 tmp->SetName(fieldname);
942
943 tmp->SetNumberOfComponents(dest_dim);
944
945 if (location == FVM_WRITER_PER_NODE) {
946 tmp->SetNumberOfTuples(f->GetNumberOfPoints());
947 // f->GetPointData()->AllocateArrays(f->GetNumberOfPoints());
948 vtkDataSetAttributes::SafeDownCast
949 (f->GetPointData())->AddArray(vtkAbstractArray::SafeDownCast(tmp));
950 }
951 else {
952 tmp->SetNumberOfTuples(f->GetNumberOfCells());
953 // f->GetCellData()->AllocateArrays(f->GetNumberOfCells());
954 /* force to AddArray, due to protection of SetArray */
955 vtkDataSetAttributes::SafeDownCast
956 (f->GetCellData())->AddArray(vtkAbstractArray::SafeDownCast(tmp));
957 }
958
959 }
960
961 BFT_REALLOC(writer->fields,
962 writer->n_fields+ 1,
963 fvm_catalyst_field_t *);
964
965 BFT_MALLOC(writer->fields[f_id], 1, fvm_catalyst_field_t);
966
967 writer->fields[f_id]->mesh_id = mesh_id;
968
969 writer->fields[f_id]->dim = dim;
970 writer->fields[f_id]->f = f;
971
972 writer->n_fields++;
973
974 return f_id;
975 }
976
977 /*----------------------------------------------------------------------------
978 * Write vertex coordinates to VTK.
979 *
980 * parameters:
981 * mesh <-- pointer to nodal mesh structure
982 * vtk_mesh <-- pointer to VTK Mesh object
983 *----------------------------------------------------------------------------*/
984
985 static void
_export_vertex_coords(const fvm_nodal_t * mesh,vtkUnstructuredGrid * vtk_mesh)986 _export_vertex_coords(const fvm_nodal_t *mesh,
987 vtkUnstructuredGrid *vtk_mesh)
988 {
989 cs_lnum_t i, j;
990 size_t stride;
991
992 double point[3];
993
994 const double *vertex_coords = mesh->vertex_coords;
995 const cs_lnum_t n_vertices = mesh->n_vertices;
996
997 /* Vertex coordinates */
998 /*--------------------*/
999
1000 stride = (size_t)(mesh->dim);
1001
1002 vtkNew<vtkPoints> points;
1003
1004 points->Allocate(mesh->n_vertices);
1005
1006 if (mesh->parent_vertex_num != NULL) {
1007 const cs_lnum_t *parent_vertex_num = mesh->parent_vertex_num;
1008 for (i = 0; i < n_vertices; i++) {
1009 for (j = 0; j < mesh->dim; j++)
1010 point[j] = vertex_coords[(parent_vertex_num[i]-1)*stride + j];
1011 for (; j < 3; j++)
1012 point[j] = 0.;
1013 points->InsertNextPoint(point[0], point[1], point[2]);
1014 }
1015 }
1016 else {
1017 for (i = 0; i < n_vertices; i++) {
1018 for (j = 0; j < mesh->dim; j++)
1019 point[j] = vertex_coords[i*stride + j];
1020 for (; j < 3; j++)
1021 point[j] = 0.;
1022 points->InsertNextPoint(point[0], point[1], point[2]);
1023 }
1024 }
1025
1026 if (mesh->global_vertex_num != NULL) {
1027
1028 const cs_gnum_t *g_vtx_num
1029 = fvm_io_num_get_global_num(mesh->global_vertex_num);
1030
1031 vtkNew<vtkIdTypeArray> g_vtx_id;
1032
1033 g_vtx_id->SetNumberOfComponents(1);
1034 g_vtx_id->SetName("GlobalNodeIds");
1035 g_vtx_id->SetNumberOfTuples(n_vertices);
1036
1037 for (i = 0; i < n_vertices; i++) {
1038 vtkIdType ii = g_vtx_num[i]-1;
1039 #if CS_PV_VERSION < 51
1040 g_vtx_id->SetTupleValue(i, &ii);
1041 #else
1042 g_vtx_id->SetTypedTuple(i, &ii);
1043 #endif
1044 }
1045
1046 vtk_mesh->GetPointData()->SetGlobalIds(g_vtx_id);
1047
1048 }
1049
1050 vtk_mesh->SetPoints(points);
1051 }
1052
1053 /*----------------------------------------------------------------------------
1054 * Write strided connectivity block to VTK.
1055 *
1056 * The connectivity on input may use 1 to n numbering, so it is shifted
1057 * by -1 here.
1058 *
1059 * parameters:
1060 * type <-- FVM element type
1061 * n_elts <-- number of elements in block
1062 * connect <-- connectivity array
1063 * vtk_mesh <-> pointer to VTK mesh object
1064 *----------------------------------------------------------------------------*/
1065
1066 static void
_write_connect_block(fvm_element_t type,cs_lnum_t n_elts,const cs_lnum_t connect[],vtkUnstructuredGrid * vtk_mesh)1067 _write_connect_block(fvm_element_t type,
1068 cs_lnum_t n_elts,
1069 const cs_lnum_t connect[],
1070 vtkUnstructuredGrid *vtk_mesh)
1071 {
1072 int vertex_order[8];
1073 cs_lnum_t i;
1074 int j;
1075
1076 vtkIdType *vtx_ids = new vtkIdType[8];
1077
1078 const int stride = fvm_nodal_n_vertices_element[type];
1079 VTKCellType vtk_type = _get_norm_elt_type(type);
1080
1081 _get_vertex_order(vtk_type, vertex_order);
1082
1083 assert(vtk_mesh != NULL);
1084
1085 for (i = 0; i < n_elts; i++) {
1086 for (j = 0; j < stride; j++)
1087 vtx_ids[j] = connect[i*stride + vertex_order[j]] - 1;
1088 vtk_mesh->InsertNextCell(vtk_type, stride, vtx_ids);
1089 }
1090
1091 delete [] vtx_ids;
1092 }
1093
1094 /*----------------------------------------------------------------------------
1095 * Write polygons from a nodal mesh to VTK.
1096 *
1097 * parameters:
1098 * export_section <-- pointer to Catalyst section helper structure
1099 * vtk_mesh <-> pointer to VTK mesh object
1100 *----------------------------------------------------------------------------*/
1101
1102 static void
_export_nodal_polygons(const fvm_nodal_section_t * section,vtkUnstructuredGrid * vtk_mesh)1103 _export_nodal_polygons(const fvm_nodal_section_t *section,
1104 vtkUnstructuredGrid *vtk_mesh)
1105 {
1106 cs_lnum_t i, j;
1107
1108 int vtx_ids_size = 8;
1109 vtkIdType *vtx_ids = NULL;
1110
1111 BFT_MALLOC(vtx_ids, vtx_ids_size, vtkIdType);
1112
1113 /* Loop on all polygonal faces */
1114 /*-----------------------------*/
1115
1116 for (i = 0; i < section->n_elements; i++) {
1117
1118 int k = 0;
1119
1120 int face_size = section->vertex_index[i+1] - section->vertex_index[i];
1121 while (vtx_ids_size < face_size) {
1122 vtx_ids_size *= 2;
1123 BFT_REALLOC(vtx_ids, vtx_ids_size, vtkIdType);
1124 }
1125
1126 for (j = section->vertex_index[i];
1127 j < section->vertex_index[i+1];
1128 j++)
1129 vtx_ids[k++] = section->vertex_num[j] - 1;
1130
1131 vtk_mesh->InsertNextCell(VTK_POLYGON, k, vtx_ids);
1132
1133 } /* End of loop on polygonal faces */
1134
1135 BFT_FREE(vtx_ids);
1136 }
1137
1138 /*----------------------------------------------------------------------------
1139 * Write polyhedra from a nodal mesh to VTK.
1140 *
1141 * parameters:
1142 * n_vertices <-- number of vertices in FVM mesh
1143 * section <-- pointer to Catalyst section helper structure
1144 * vtk_mesh <-> VTK mesh object
1145 *----------------------------------------------------------------------------*/
1146
1147 static void
_export_nodal_polyhedra(cs_lnum_t n_vertices,const fvm_nodal_section_t * section,vtkUnstructuredGrid * vtk_mesh)1148 _export_nodal_polyhedra(cs_lnum_t n_vertices,
1149 const fvm_nodal_section_t *section,
1150 vtkUnstructuredGrid *vtk_mesh)
1151 {
1152 int face_sgn;
1153 int buf_size = 8;
1154 cs_lnum_t i, j, k, l;
1155
1156 cs_lnum_t face_length, face_id;
1157
1158 int *vtx_marker = NULL;
1159 vtkIdType *face_array = NULL;
1160 vtkIdType *vtx_ids = NULL;
1161
1162 BFT_MALLOC(vtx_marker, n_vertices, int);
1163 BFT_MALLOC(face_array, buf_size, vtkIdType);
1164 BFT_MALLOC(vtx_ids, buf_size, vtkIdType);
1165
1166 for (i = 0; i < n_vertices; i++)
1167 vtx_marker[i] = -1;
1168
1169 /* Write cell/vertex connectivity */
1170 /*--------------------------------*/
1171
1172 for (i = 0; i < section->n_elements; i++) {
1173
1174 int m = 0;
1175 int cell_vtx_count = 0;
1176 int n_elt_faces = section->face_index[i+1] - section->face_index[i];
1177
1178 for (j = section->face_index[i];
1179 j < section->face_index[i+1];
1180 j++) {
1181
1182 if (section->face_num[j] > 0) {
1183 face_id = section->face_num[j] - 1;
1184 face_sgn = 1;
1185 }
1186 else {
1187 face_id = -section->face_num[j] - 1;
1188 face_sgn = -1;
1189 }
1190
1191 face_length = ( section->vertex_index[face_id+1]
1192 - section->vertex_index[face_id]);
1193
1194 while (m + face_length + 1 > buf_size) {
1195 buf_size *= 2;
1196 BFT_REALLOC(face_array, buf_size, vtkIdType);
1197 BFT_REALLOC(vtx_ids, buf_size, vtkIdType);
1198 }
1199
1200 face_array[m++] = face_length;
1201
1202 for (k = 0; k < face_length; k++) {
1203 int vtx_id;
1204 l = section->vertex_index[face_id]
1205 + (face_length + (k*face_sgn))%face_length;
1206 vtx_id = section->vertex_num[l] - 1;
1207 face_array[m++] = vtx_id;
1208 if (vtx_marker[vtx_id] < i) {
1209 vtx_ids[cell_vtx_count] = vtx_id;
1210 vtx_marker[vtx_id] = i;
1211 cell_vtx_count += 1;
1212 }
1213 }
1214
1215 } /* End of loop on cell faces */
1216
1217 vtk_mesh->InsertNextCell(VTK_POLYHEDRON, cell_vtx_count,
1218 vtx_ids, n_elt_faces, face_array);
1219
1220 } /* End of loop on polyhedral cells */
1221
1222 BFT_FREE(vtx_ids);
1223 BFT_FREE(face_array);
1224 BFT_FREE(vtx_marker);
1225 }
1226
1227 /*----------------------------------------------------------------------------
1228 * Write field values associated with nodal values of a nodal mesh to VTK.
1229 *
1230 * Output fields ar either scalar or 3d vectors or scalars, and are
1231 * non interlaced. Input arrays may be less than 2d, in which case the z
1232 * values are set to 0, and may be interlaced or not.
1233 *
1234 * parameters:
1235 * mesh <-- pointer to nodal mesh structure
1236 * dim <-- field dimension
1237 * interlace <-- indicates if field in memory is interlaced
1238 * n_parent_lists <-- indicates if field values are to be obtained
1239 * directly through the local entity index (when 0) or
1240 * through the parent entity numbers (when 1 or more)
1241 * parent_num_shift <-- parent list to common number index shifts;
1242 * size: n_parent_lists
1243 * datatype <-- input data type (output is real)
1244 * field_values <-- array of associated field value arrays
1245 * f <-- associated vtkUnstructuredGrid object
1246 *----------------------------------------------------------------------------*/
1247
1248 static void
_export_field_values_n(const fvm_nodal_t * mesh,const char * fieldname,int dim,cs_interlace_t interlace,int n_parent_lists,const cs_lnum_t parent_num_shift[],cs_datatype_t datatype,const void * const field_values[],vtkUnstructuredGrid * f)1249 _export_field_values_n(const fvm_nodal_t *mesh,
1250 const char *fieldname,
1251 int dim,
1252 cs_interlace_t interlace,
1253 int n_parent_lists,
1254 const cs_lnum_t parent_num_shift[],
1255 cs_datatype_t datatype,
1256 const void *const field_values[],
1257 vtkUnstructuredGrid *f)
1258 {
1259 assert(f != NULL);
1260
1261 double *values = NULL;
1262
1263 const int dest_dim = (dim == 6) ? 9 : dim;
1264
1265 values = vtkDoubleArray::SafeDownCast
1266 (f->GetPointData()->GetArray(fieldname))
1267 ->WritePointer(0, dest_dim*f->GetNumberOfPoints());
1268
1269 fvm_convert_array(dim,
1270 0,
1271 dest_dim,
1272 0,
1273 mesh->n_vertices,
1274 interlace,
1275 datatype,
1276 CS_DOUBLE,
1277 n_parent_lists,
1278 parent_num_shift,
1279 mesh->parent_vertex_num,
1280 field_values,
1281 values);
1282
1283 /* Special case for symmetric tensors */
1284
1285 if (dim == 6) {
1286 for (cs_lnum_t i = 0; i < mesh->n_vertices; i++) {
1287 values[9*i + 8] = values[9*i + 2];
1288 values[9*i + 7] = values[9*i + 4];
1289 values[9*i + 6] = values[9*i + 5];
1290 values[9*i + 4] = values[9*i + 1];
1291 values[9*i + 2] = values[9*i + 5];
1292 values[9*i + 1] = values[9*i + 3];
1293 values[9*i + 5] = values[9*i + 7];
1294 }
1295 }
1296
1297 }
1298
1299 /*----------------------------------------------------------------------------
1300 * Write field values associated with element values of a nodal mesh to VTK.
1301 *
1302 * Output fields are non interlaced. Input arrays may be interlaced or not.
1303 *
1304 * parameters:
1305 * mesh <-- pointer to nodal mesh structure
1306 * dim <-- field dimension
1307 * interlace <-- indicates if field in memory is interlaced
1308 * n_parent_lists <-- indicates if field values are to be obtained
1309 * directly through the local entity index (when 0) or
1310 * through the parent entity numbers (when 1 or more)
1311 * parent_num_shift <-- parent list to common number index shifts;
1312 * size: n_parent_lists
1313 * datatype <-- indicates the data type of (source) field values
1314 * field_values <-- array of associated field value arrays
1315 * f <-- associated object handle
1316 *----------------------------------------------------------------------------*/
1317
1318 static void
_export_field_values_e(const fvm_nodal_t * mesh,const char * fieldname,int dim,cs_interlace_t interlace,int n_parent_lists,const cs_lnum_t parent_num_shift[],cs_datatype_t datatype,const void * const field_values[],vtkUnstructuredGrid * f)1319 _export_field_values_e(const fvm_nodal_t *mesh,
1320 const char *fieldname,
1321 int dim,
1322 cs_interlace_t interlace,
1323 int n_parent_lists,
1324 const cs_lnum_t parent_num_shift[],
1325 cs_datatype_t datatype,
1326 const void *const field_values[],
1327 vtkUnstructuredGrid *f)
1328 {
1329 assert(f != NULL);
1330
1331 int section_id;
1332 double *values = NULL;
1333
1334 const int dest_dim = (dim == 6) ? 9 : dim;
1335
1336 values = vtkDoubleArray::SafeDownCast
1337 (f->GetCellData()->GetArray(fieldname))
1338 ->WritePointer(0, dest_dim*f->GetNumberOfCells());
1339
1340 /* Distribute partition to block values */
1341
1342 cs_lnum_t start_id = 0;
1343 cs_lnum_t src_shift = 0;
1344
1345 /* loop on sections which should be appended */
1346
1347 const int elt_dim = fvm_nodal_get_max_entity_dim(mesh);
1348
1349 for (section_id = 0; section_id < mesh->n_sections; section_id++) {
1350
1351 const fvm_nodal_section_t *section = mesh->sections[section_id];
1352
1353 if (section->entity_dim < elt_dim)
1354 continue;
1355
1356 assert(values != NULL || section->n_elements == 0);
1357
1358 fvm_convert_array(dim,
1359 0,
1360 dest_dim,
1361 src_shift,
1362 section->n_elements + src_shift,
1363 interlace,
1364 datatype,
1365 CS_DOUBLE,
1366 n_parent_lists,
1367 parent_num_shift,
1368 section->parent_element_num,
1369 field_values,
1370 values + start_id);
1371
1372 start_id += section->n_elements*dest_dim;
1373 if (n_parent_lists == 0)
1374 src_shift += section->n_elements;
1375
1376 }
1377
1378 /* Special case for symmetric tensors */
1379
1380 if (dim == 6) {
1381
1382 cs_lnum_t n_elts = f->GetNumberOfCells();
1383
1384 assert(values != NULL || n_elts == 0);
1385
1386 for (cs_lnum_t i = 0; i < n_elts; i++) {
1387 values[9*i + 8] = values[9*i + 2];
1388 values[9*i + 7] = values[9*i + 4];
1389 values[9*i + 6] = values[9*i + 5];
1390 values[9*i + 4] = values[9*i + 1];
1391 values[9*i + 2] = values[9*i + 5];
1392 values[9*i + 1] = values[9*i + 3];
1393 values[9*i + 5] = values[9*i + 7];
1394 }
1395 }
1396 }
1397
1398 /*============================================================================
1399 * Public function definitions
1400 *============================================================================*/
1401
1402 BEGIN_C_DECLS
1403
1404 /*----------------------------------------------------------------------------
1405 * Initialize FVM to Catalyst object writer.
1406 *
1407 * Options are:
1408 * private_comm use private MPI communicator (default: false)
1409 * names=<fmt> use same naming rules as <fmt> format
1410 * (default: ensight)
1411 * input_name=<name> define input name (default: writer name)
1412 *
1413 * parameters:
1414 * name <-- base output case name.
1415 * options <-- whitespace separaed, lowercase options list
1416 * time_dependecy <-- indicates if and how meshes will change with time
1417 * comm <-- associated MPI communicator.
1418 *
1419 * returns:
1420 * pointer to opaque Catalyst writer structure.
1421 *----------------------------------------------------------------------------*/
1422
1423 #if defined(HAVE_MPI)
1424 void *
fvm_to_catalyst_init_writer(const char * name,const char * path,const char * options,fvm_writer_time_dep_t time_dependency,MPI_Comm comm)1425 fvm_to_catalyst_init_writer(const char *name,
1426 const char *path,
1427 const char *options,
1428 fvm_writer_time_dep_t time_dependency,
1429 MPI_Comm comm)
1430 #else
1431 void *
1432 fvm_to_catalyst_init_writer(const char *name,
1433 const char *path,
1434 const char *options,
1435 fvm_writer_time_dep_t time_dependency)
1436 #endif
1437 {
1438 CS_UNUSED(path);
1439
1440 fvm_to_catalyst_t *w = NULL;
1441
1442 bool private_comm = false;
1443
1444 /* Initialize writer */
1445
1446 BFT_MALLOC(w, 1, fvm_to_catalyst_t);
1447
1448 w->rank = 0;
1449 w->n_ranks = 1;
1450
1451 w->time_dependency = time_dependency;
1452
1453 w->mb = vtkMultiBlockDataSet::New();
1454
1455 w->n_fields = 0;
1456 w->fields = NULL;
1457
1458 w->time_step = -1;
1459 w->time_value = 0.0;
1460
1461 w->ensight_names = true;
1462 w->input_name = NULL;
1463
1464 /* Writer name */
1465
1466 if (name != NULL) {
1467 BFT_MALLOC(w->name, strlen(name) + 1, char);
1468 strcpy(w->name, name);
1469 }
1470 else {
1471 const char _name[] = "catalyst";
1472 BFT_MALLOC(w->name, strlen(_name) + 1, char);
1473 strcpy(w->name, _name);
1474 }
1475
1476 /* Parse options */
1477
1478 if (options != NULL) {
1479
1480 int i1, i2, l_opt;
1481 int l_tot = strlen(options);
1482
1483 i1 = 0; i2 = 0;
1484 while (i1 < l_tot) {
1485
1486 for (i2 = i1; i2 < l_tot && options[i2] != ' '; i2++);
1487 l_opt = i2 - i1;
1488
1489 if ((l_opt == 12) && (strncmp(options + i1, "private_comm", l_opt) == 0))
1490 private_comm = true;
1491 else if ((l_opt > 6) && (strncmp(options + i1, "names=", 6) == 0)) {
1492 if ((l_opt == 6+7) && (strncmp(options + i1 + 6, "ensight", 7) == 0))
1493 w->ensight_names = true;
1494 else
1495 w->ensight_names = false;
1496 }
1497 else if ((l_opt > 11) && (strncmp(options + i1, "input_name=", 11) == 0)) {
1498 int l = strlen(options + i1 + 11);
1499 BFT_MALLOC(w->input_name, l+1, char);
1500 strncpy(w->input_name, options + i1 + 11, l);
1501 w->input_name[l] = '\0';
1502 }
1503
1504 for (i1 = i2 + 1; i1 < l_tot && options[i1] == ' '; i1++);
1505
1506 }
1507
1508 }
1509
1510 #if CS_PV_VERSION < 55
1511 if (w->input_name == NULL) {
1512 char n[] = "input";
1513 BFT_MALLOC(w->input_name, strlen(n)+1, char);
1514 strcpy(w->input_name, n);
1515 }
1516 #endif
1517
1518 /* Ensure coprocessor is built */
1519
1520 _n_writers += 1;
1521
1522 #if defined(HAVE_MPI)
1523 _init_coprocessor(private_comm, comm);
1524 #else
1525 _init_coprocessor();
1526 #endif
1527
1528 /* Parallel parameters */
1529
1530 #if defined(HAVE_MPI)
1531 if (_comm != MPI_COMM_NULL) {
1532 MPI_Comm_rank(_comm, &(w->rank));
1533 MPI_Comm_size(_comm, &(w->n_ranks));
1534 }
1535 #endif
1536
1537 if (_n_scripts < 1)
1538 _add_dir_scripts(".");
1539
1540 w->modified = true;
1541
1542 return w;
1543 }
1544
1545 /*----------------------------------------------------------------------------
1546 * Finalize FVM to Catalyst object writer.
1547 *
1548 * parameters:
1549 * this_writer_p <-- pointer to opaque writer structure.
1550 *
1551 * returns:
1552 * NULL pointer
1553 *----------------------------------------------------------------------------*/
1554
1555 void *
fvm_to_catalyst_finalize_writer(void * this_writer_p)1556 fvm_to_catalyst_finalize_writer(void *this_writer_p)
1557 {
1558 int i;
1559
1560 fvm_to_catalyst_t *w = (fvm_to_catalyst_t *)this_writer_p;
1561
1562 assert(w != NULL);
1563
1564 /* Write output if not done already */
1565
1566 fvm_to_catalyst_flush(this_writer_p);
1567
1568 /* Free structures */
1569
1570 BFT_FREE(w->name);
1571
1572 /* Free vtkUnstructuredGrid and field structures
1573 (reference counters should go to 0) */
1574
1575 _free_coprocessor();
1576
1577 _n_writers -= 1;
1578
1579 w->mb->Delete();
1580
1581 for (i = 0; i < w->n_fields; i++) {
1582 w->fields[i]->f = NULL; // delete w->fields[i]->f;
1583 BFT_FREE(w->fields[i]);
1584 }
1585
1586 BFT_FREE(w->fields);
1587
1588 /* Free fvm_to_catalyst_t structure */
1589
1590 BFT_FREE(w);
1591
1592 return NULL;
1593 }
1594
1595 /*----------------------------------------------------------------------------
1596 * Associate new time step with an Catalyst geometry.
1597 *
1598 * parameters:
1599 * this_writer_p <-- pointer to associated writer
1600 * time_step <-- time step number
1601 * time_value <-- time_value number
1602 *----------------------------------------------------------------------------*/
1603
1604 void
fvm_to_catalyst_set_mesh_time(void * this_writer_p,int time_step,double time_value)1605 fvm_to_catalyst_set_mesh_time(void *this_writer_p,
1606 int time_step,
1607 double time_value)
1608 {
1609 fvm_to_catalyst_t *w = (fvm_to_catalyst_t *)this_writer_p;
1610
1611 int _time_step = (time_step > -1) ? time_step : 0;
1612 double _time_value = (time_value > 0.0) ? time_value : 0.0;
1613
1614 if (_time_step > w->time_step) {
1615 if ( w->time_dependency == FVM_WRITER_TRANSIENT_CONNECT
1616 && _time_step > w->time_step) {
1617 for (int i = 0; i < w->n_fields; i++) {
1618 w->fields[i]->f = NULL; // delete w->fields[i]->f;
1619 BFT_FREE(w->fields[i]);
1620 }
1621 BFT_FREE(w->fields);
1622 w->n_fields = 0;
1623 w->mb->Delete();
1624 w->mb = vtkMultiBlockDataSet::New();
1625 }
1626 w->time_step = _time_step;
1627 assert(time_value >= w->time_value);
1628 w->time_value = _time_value;
1629 }
1630 }
1631
1632 /*----------------------------------------------------------------------------
1633 * Write nodal mesh to a Catalyst object
1634 *
1635 * parameters:
1636 * this_writer_p <-- pointer to associated writer
1637 * mesh <-- pointer to nodal mesh structure that should be written
1638 *----------------------------------------------------------------------------*/
1639
1640 void
fvm_to_catalyst_export_nodal(void * this_writer_p,const fvm_nodal_t * mesh)1641 fvm_to_catalyst_export_nodal(void *this_writer_p,
1642 const fvm_nodal_t *mesh)
1643 {
1644 int mesh_id, section_id;
1645
1646 fvm_to_catalyst_t *w = (fvm_to_catalyst_t *)this_writer_p;
1647
1648 const int elt_dim = fvm_nodal_get_max_entity_dim(mesh);
1649
1650 /* Initialization */
1651 /*----------------*/
1652
1653 /* Get matching mesh */
1654
1655 mesh_id = _get_catalyst_mesh_id(w, mesh->name);
1656
1657 if (mesh_id < 0)
1658 mesh_id = _add_catalyst_mesh(w, mesh);
1659
1660 vtkUnstructuredGrid *ugrid
1661 = vtkUnstructuredGrid::SafeDownCast
1662 (vtkDataSet::SafeDownCast(w->mb->GetBlock(mesh_id)));
1663
1664 /* Vertex coordinates */
1665 /*--------------------*/
1666
1667 _export_vertex_coords(mesh, ugrid);
1668
1669 /* Element connectivity size */
1670 /*---------------------------*/
1671
1672 cs_lnum_t n_elts = 0;
1673
1674 for (section_id = 0; section_id < mesh->n_sections; section_id++) {
1675
1676 const fvm_nodal_section_t *section = mesh->sections[section_id];
1677
1678 if (section->entity_dim < elt_dim)
1679 continue;
1680
1681 n_elts += section->n_elements;
1682
1683 } /* End of loop on sections */
1684
1685 if (n_elts > 0)
1686 ugrid->Allocate(n_elts);
1687
1688 /* Element connectivity */
1689 /*----------------------*/
1690
1691 for (section_id = 0; section_id < mesh->n_sections; section_id++) {
1692
1693 const fvm_nodal_section_t *section = mesh->sections[section_id];
1694
1695 if (section->entity_dim < elt_dim)
1696 continue;
1697
1698 if (section->stride > 0)
1699 _write_connect_block(section->type,
1700 section->n_elements,
1701 section->vertex_num,
1702 ugrid);
1703
1704 else if (section->type == FVM_FACE_POLY)
1705 _export_nodal_polygons(section, ugrid);
1706
1707 else if (section->type == FVM_CELL_POLY)
1708 _export_nodal_polyhedra(mesh->n_vertices, section, ugrid);
1709
1710 } /* End of loop on sections */
1711
1712 w->modified = true;
1713 }
1714
1715 /*----------------------------------------------------------------------------
1716 * Write field associated with a nodal mesh to a Catalyst object.
1717 *
1718 * Assigning a negative value to the time step indicates a time-independent
1719 * field (in which case the time_value argument is unused).
1720 *
1721 * parameters:
1722 * this_writer_p <-- pointer to associated writer
1723 * mesh <-- pointer to associated nodal mesh structure
1724 * name <-- variable name
1725 * location <-- variable definition location (nodes or elements)
1726 * dimension <-- variable dimension (0: constant, 1: scalar,
1727 * 3: vector, 6: sym. tensor, 9: asym. tensor)
1728 * interlace <-- indicates if variable in memory is interlaced
1729 * n_parent_lists <-- indicates if variable values are to be obtained
1730 * directly through the local entity index (when 0) or
1731 * through the parent entity numbers (when 1 or more)
1732 * parent_num_shift <-- parent number to value array index shifts;
1733 * size: n_parent_lists
1734 * datatype <-- indicates the data type of (source) field values
1735 * time_step <-- number of the current time step
1736 * time_value <-- associated time value
1737 * field_values <-- array of associated field value arrays
1738 *----------------------------------------------------------------------------*/
1739
1740 void
fvm_to_catalyst_export_field(void * this_writer_p,const fvm_nodal_t * mesh,const char * name,fvm_writer_var_loc_t location,int dimension,cs_interlace_t interlace,int n_parent_lists,const cs_lnum_t parent_num_shift[],cs_datatype_t datatype,int time_step,double time_value,const void * const field_values[])1741 fvm_to_catalyst_export_field(void *this_writer_p,
1742 const fvm_nodal_t *mesh,
1743 const char *name,
1744 fvm_writer_var_loc_t location,
1745 int dimension,
1746 cs_interlace_t interlace,
1747 int n_parent_lists,
1748 const cs_lnum_t parent_num_shift[],
1749 cs_datatype_t datatype,
1750 int time_step,
1751 double time_value,
1752 const void *const field_values[])
1753 {
1754 int mesh_id, field_id;
1755 char _name[128];
1756
1757 fvm_to_catalyst_t *w = (fvm_to_catalyst_t *)this_writer_p;
1758
1759 /* Initialization */
1760 /*----------------*/
1761
1762 mesh_id = _get_catalyst_mesh_id(w, mesh->name);
1763
1764 strncpy(_name, name, 127);
1765 _name[127] = '\0';
1766 if (w->ensight_names) {
1767 for (int i = 0; i < 127 && _name[i] != '\0'; i++) {
1768 switch (_name[i]) {
1769 case '(':
1770 case ')':
1771 case ']':
1772 case '[':
1773 case '+':
1774 case '-':
1775 case '@':
1776 case ' ':
1777 case '\t':
1778 case '!':
1779 case '#':
1780 case '*':
1781 case '^':
1782 case '$':
1783 case '/':
1784 _name[i] = '_';
1785 break;
1786 default:
1787 break;
1788 }
1789 if (_name[i] == ' ')
1790 _name[i] = '_';
1791 }
1792 }
1793
1794 if (mesh_id < 0) {
1795 mesh_id = _add_catalyst_mesh(w, mesh);
1796 fvm_to_catalyst_export_nodal(w, mesh);
1797 }
1798
1799 int _time_step = (time_step > -1) ? time_step : 0;
1800 double _time_value = (time_value > 0.0) ? time_value : 0.0;
1801 if (_time_step > w->time_step) {
1802 w->time_step = _time_step;
1803 assert(time_value >= w->time_value);
1804 w->time_value = _time_value;
1805 }
1806
1807 /* Get field id */
1808
1809 field_id = _get_catalyst_field_id(w,
1810 _name,
1811 mesh_id,
1812 dimension,
1813 datatype,
1814 location);
1815
1816 if (field_id < 0)
1817 field_id = _add_catalyst_field(w,
1818 _name,
1819 mesh_id,
1820 dimension,
1821 datatype,
1822 location);
1823
1824 vtkUnstructuredGrid *f = w->fields[field_id]->f;
1825
1826 /* Per node variable */
1827 /*-------------------*/
1828
1829 if (location == FVM_WRITER_PER_NODE)
1830 _export_field_values_n(mesh,
1831 _name,
1832 dimension,
1833 interlace,
1834 n_parent_lists,
1835 parent_num_shift,
1836 datatype,
1837 field_values,
1838 f);
1839
1840
1841 /* Per element variable */
1842 /*----------------------*/
1843
1844 else if (location == FVM_WRITER_PER_ELEMENT)
1845 _export_field_values_e(mesh,
1846 _name,
1847 dimension,
1848 interlace,
1849 n_parent_lists,
1850 parent_num_shift,
1851 datatype,
1852 field_values,
1853 f);
1854
1855 /* Update field status */
1856 /*---------------------*/
1857
1858 fvm_to_catalyst_set_mesh_time(w, time_step, time_value);
1859
1860 w->modified = true;
1861 }
1862
1863 /*----------------------------------------------------------------------------
1864 * Flush files associated with a given writer.
1865 *
1866 * In this case, the effective call to coprocessing is done.
1867 *
1868 * parameters:
1869 * this_writer_p <-- pointer to associated writer
1870 *----------------------------------------------------------------------------*/
1871
1872 void
fvm_to_catalyst_flush(void * this_writer_p)1873 fvm_to_catalyst_flush(void *this_writer_p)
1874 {
1875 fvm_to_catalyst_t *w = (fvm_to_catalyst_t *)this_writer_p;
1876
1877 vtkNew<vtkCPDataDescription> dataDescription;
1878 if (w->input_name != NULL)
1879 dataDescription->AddInput(w->input_name);
1880 else
1881 dataDescription->AddInput(w->name);
1882 dataDescription->SetTimeData(w->time_value, w->time_step);
1883
1884 if (_processor->RequestDataDescription(dataDescription) != 0 && w->modified) {
1885 int n = dataDescription->GetNumberOfInputDescriptions();
1886 if (n == 1)
1887 dataDescription->GetInputDescription(0)->SetGrid(w->mb);
1888 else {
1889 if (w->input_name != NULL)
1890 dataDescription->GetInputDescriptionByName(w->input_name)->SetGrid(w->mb);
1891 else
1892 dataDescription->GetInputDescriptionByName(w->name)->SetGrid(w->mb);
1893 }
1894
1895 _processor->CoProcess(dataDescription);
1896 w->modified = false;
1897 }
1898 }
1899
1900 /*----------------------------------------------------------------------------*/
1901
1902 END_C_DECLS
1903
1904 #endif /* HAVE_CATALYST */
1905