SFEMaNS  version 5.3
Reference documentation for SFEMaNS
plot_vtk.f90
Go to the documentation of this file.
1 MODULE vtk_viz
2 #include "petsc/finclude/petsc.h"
3  USE petsc
4  PUBLIC :: create_pvd_file, check_list,&
5  make_vtu_file_arpack
6  PUBLIC :: create_vtu_file_axi3d
7  PUBLIC :: make_vtu_file_3d, create_xml_vtu_scal_file, create_xml_vtu_vect_file
8  !PUBLIC :: make_vtu_file_scalar_2D
9  !PUBLIC :: make_vtu_file_axi3D, create_vtk_file
10  PRIVATE
11 CONTAINS
12 
13  SUBROUTINE check_list(communicator, file_list, check)
14  IMPLICIT NONE
15  CHARACTER(LEN=200), DIMENSION(:), POINTER :: file_list
16  CHARACTER(LEN=200), DIMENSION(:), POINTER :: dummy_list
17  INTEGER, DIMENSION(SIZE(file_list)) :: check_mylist
18  INTEGER :: check, n, count
19 !#include "petsc/finclude/petsc.h"
20  mpi_comm :: communicator
21  petscmpiint :: rank, nb_procs
22  petscerrorcode :: ierr
23  CALL mpi_comm_rank(communicator, rank, ierr)
24  CALL mpi_comm_size(communicator, nb_procs, ierr)
25 
26  CALL mpi_allgather(check, 1, mpi_integer, check_mylist, 1, &
27  mpi_integer, communicator, ierr)
28 
29  count = 0
30  DO n = 1, SIZE(file_list)
31  IF (check_mylist(n)==0) cycle
32  count = count + 1
33  END DO
34  ALLOCATE(dummy_list(count))
35  count = 0
36  DO n = 1, SIZE(file_list)
37  IF (check_mylist(n)==0) cycle
38  count = count + 1
39  dummy_list(count) = file_list(n)
40  END DO
41  DEALLOCATE(file_list)
42  ALLOCATE(file_list(count))
43  file_list = dummy_list
44  END SUBROUTINE check_list
45 
46  SUBROUTINE create_pvd_file(file_list, file_header, time_step, what)
47  IMPLICIT NONE
48  CHARACTER(*), DIMENSION(:), INTENT(IN) :: file_list
49  CHARACTER(*), INTENT(IN) :: file_header, what
50  INTEGER, INTENT(IN) :: time_step
51  INTEGER :: unit_file=789, j
52  CHARACTER(len=5) :: tit, tit_part
53 
54  IF (what=='new') THEN
55  OPEN (unit=unit_file, file=file_header//'.pvd', form = 'formatted', &
56  access = 'append', status = 'replace')
57  WRITE(unit_file, '(A)') '<?xml version="1.0"?>'
58  WRITE(unit_file, '(A)') '<VTKFile type="Collection" version="0.1"'// &
59  ' byte_order="LittleEndian" compressor="vtkZLibDataCompressor">'
60  WRITE(unit_file, '(A)') '<Collection>'
61  ELSE
62  OPEN (unit=unit_file, file=file_header//'.pvd', form = 'formatted', &
63  access = 'append', status = 'old')
64  backspace(unit_file)
65  backspace(unit_file)
66  END IF
67  WRITE(tit,'(I5)') time_step
68  DO j = 1, SIZE(file_list)
69  WRITE(tit_part,'(I5)') j
70  WRITE(unit_file,'(A)') '<DataSet timestep="'//trim(adjustl(tit))//&
71  '" group="" part="'// &
72  trim(adjustl(tit_part))//'" file="./'//trim(adjustl(file_list(j)))//&
73  '.vtu'//'"/>'
74  END DO
75  WRITE(unit_file, '(A)') '</Collection>'
76  WRITE(unit_file, '(A)') '</VTKFile>'
77 
78  CLOSE(unit_file) !Ecriture pour paraview
79  END SUBROUTINE create_pvd_file
80 
81  SUBROUTINE create_xml_vtu_scal_file(field, mesh, file_name, field_name)
82  USE def_type_mesh
83  USE input_data
84  USE zlib_base64
85  IMPLICIT NONE
86  TYPE(mesh_type) :: mesh
87  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: field
88  CHARACTER(*), INTENT(IN) :: file_name, field_name
89  INTEGER :: unit_file=789, m, n, type_cell
90  REAL(KIND=4), DIMENSION(3*mesh%np) :: r4_threed_xml_field
91  INTEGER(KIND=4), DIMENSION(mesh%gauss%n_w*mesh%me) :: i4_threed_xml_field
92  INTEGER(KIND=4), DIMENSION(mesh%me) :: i4_xml_field
93  INTEGER(KIND=1), DIMENSION(mesh%me) :: i1_xml_field
94  CHARACTER(LEN=200) :: ascii_or_binary
95 
96  IF (SIZE(field)==0) RETURN
97 
98  IF (inputs%if_xml) THEN
99  ascii_or_binary = 'binary'
100  OPEN (unit=unit_file, file=file_name//'.vtu', status = 'unknown')
101  WRITE(unit_file,'(A)') '<?xml version="1.0" ?>'
102  WRITE(unit_file,'(A)', advance="no") '<VTKFile type="UnstructuredGrid" version="0.1" '
103  WRITE(unit_file,'(A)') 'compressor="vtkZLibDataCompressor" byte_order="LittleEndian">'
104  ELSE
105  ascii_or_binary = 'ascii'
106  OPEN (unit=unit_file, file=file_name//'.vtu',&
107  form = 'formatted', status = 'unknown')
108  WRITE(unit_file,'(A)') '<VTKFile type="UnstructuredGrid" version="0.1"'// &
109  ' byte_order="LittleEndian">'
110  END IF
111 
112  WRITE(unit_file,'(A)') '<UnstructuredGrid>'
113  WRITE(unit_file,'(A,I9,A,I9,A)') '<Piece NumberOfPoints="', mesh%np, &
114  '" NumberOfCells="', mesh%me, '">'
115 
116  !===PointData Block
117  WRITE(unit_file,'(A)') '<PointData Scalars="truc">'
118  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="'&
119  //trim(adjustl(field_name))//'" format="'//trim(adjustl(ascii_or_binary))//'">'
120  IF (inputs%if_xml) THEN
121  CALL write_compressed_block(unit_file, REAL(field(1:mesh%np),4))
122  ELSE
123  DO n = 1, mesh%np
124  WRITE(unit_file,'(e14.7)') field(n)
125  ENDDO
126  END IF
127  WRITE(unit_file,'(A)') '</DataArray>'
128  WRITE(unit_file,'(A)') '</PointData>'
129  !===End of PointData Block
130 
131  !===CellData Block
132  WRITE(unit_file,'(A)') '<CellData>'
133  WRITE(unit_file,'(A)') '</CellData>'
134  !===End of CellData Block
135 
136  !===Points Block
137  WRITE(unit_file,'(A)') '<Points>'
138  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="Points" '//&
139  'NumberOfComponents="3" format="'//trim(adjustl(ascii_or_binary))//'">'
140  IF (inputs%if_xml) THEN
141  DO n = 1, mesh%np
142  r4_threed_xml_field(3*(n-1)+1) = REAL(mesh%rr(1,n),4)
143  r4_threed_xml_field(3*(n-1)+2) = REAL(0.,4)
144  r4_threed_xml_field(3*(n-1)+3) = REAL(mesh%rr(2,n),4)
145  END DO
146  CALL write_compressed_block(unit_file, r4_threed_xml_field(1:3*mesh%np))
147  ELSE
148  DO n = 1, mesh%np
149  WRITE(unit_file,'(3(e14.7,x))') mesh%rr(1,n), 0.d0 , mesh%rr(3,n)
150  END DO
151  END IF
152  WRITE(unit_file,'(A)') '</DataArray>'
153  WRITE(unit_file,'(A)') '</Points>'
154  !===End of Points Block
155 
156  !===Cells Block
157  IF (mesh%gauss%n_w==3) THEN
158  type_cell = 5
159  ELSE IF (mesh%gauss%n_w==6) THEN
160  type_cell = 22
161  END IF
162  WRITE(unit_file,'(A)') '<Cells>'
163  WRITE(unit_file,'(A)') '<DataArray type="Int32" Name="connectivity" format="'&
164  //trim(adjustl(ascii_or_binary))//'">'
165  IF (inputs%if_xml) THEN
166  DO m = 1, mesh%me
167  i4_threed_xml_field(mesh%gauss%n_w*(m-1)+1) = int(mesh%jj(1,m)-1,4)
168  i4_threed_xml_field(mesh%gauss%n_w*(m-1)+2) = int(mesh%jj(2,m)-1,4)
169  i4_threed_xml_field(mesh%gauss%n_w*(m-1)+3) = int(mesh%jj(3,m)-1,4)
170  IF (type_cell==22) THEN
171  i4_threed_xml_field(mesh%gauss%n_w*(m-1)+4) = int(mesh%jj(6,m)-1,4)
172  i4_threed_xml_field(mesh%gauss%n_w*(m-1)+5) = int(mesh%jj(4,m)-1,4)
173  i4_threed_xml_field(mesh%gauss%n_w*(m-1)+6) = int(mesh%jj(5,m)-1,4)
174  END IF
175  END DO
176  CALL write_compressed_block(unit_file, i4_threed_xml_field)
177  ELSE
178  DO m = 1, mesh%me
179  WRITE(unit_file,'(3(I8,1x))') mesh%jj(1:3,m)-1
180  IF (type_cell==22) THEN
181  WRITE(unit_file,'(3(I8,1x))') mesh%jj(6,m)-1 , mesh%jj(4,m)-1 , mesh%jj(5,m)-1
182  END IF
183  END DO
184  END IF
185  WRITE(unit_file,'(A)') '</DataArray>'
186 
187  WRITE(unit_file,'(A)') '<DataArray type="Int32" Name="offsets" format="'&
188  //trim(adjustl(ascii_or_binary))//'">'
189  IF (inputs%if_xml) THEN
190  DO m = 1, mesh%me
191  i4_xml_field(m) = int(mesh%gauss%n_w*m,4)
192  END DO
193  CALL write_compressed_block(unit_file, i4_xml_field)
194  ELSE
195  DO m = 1, mesh%me
196  WRITE(unit_file,'(I8)') m*mesh%gauss%n_w
197  END DO
198  END IF
199  WRITE(unit_file,'(A)') '</DataArray>'
200 
201  WRITE(unit_file,'(A)') '<DataArray type="UInt8" Name="types" format="'&
202  //trim(adjustl(ascii_or_binary))//'">'
203  IF (inputs%if_xml) THEN
204  DO m = 1, mesh%me
205  i1_xml_field(m) = int(type_cell,1)
206  END DO
207  CALL write_compressed_block(unit_file, i1_xml_field)
208  ELSE
209  DO m = 1, mesh%me
210  WRITE(unit_file,'(I8)') type_cell
211  END DO
212  END IF
213  WRITE(unit_file,'(A)') '</DataArray>'
214  WRITE(unit_file,'(A)') '</Cells>'
215  !===End of Cells Block
216 
217  WRITE(unit_file,'(A)') '</Piece>'
218  WRITE(unit_file,'(A)') '</UnstructuredGrid>'
219  WRITE(unit_file,'(A)') '</VTKFile>'
220 
221  CLOSE(unit_file)
222  END SUBROUTINE create_xml_vtu_scal_file
223 
224  SUBROUTINE create_xml_vtu_vect_file(field, mesh, file_name, opt_st)
225  USE def_type_mesh
226  USE input_data
227  USE zlib_base64
228  IMPLICIT NONE
229  TYPE(mesh_type) :: mesh
230  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: field
231  CHARACTER(*), INTENT(IN) :: file_name
232  INTEGER :: unit_file=789, m, n, type_cell
233  CHARACTER(*), OPTIONAL, INTENT(IN) :: opt_st
234  CHARACTER(LEN=200) :: field_name
235  REAL(KIND=4), DIMENSION(3*mesh%np) :: r4_threed_xml_field
236  INTEGER(KIND=4), DIMENSION(mesh%gauss%n_w*mesh%me) :: i4_threed_xml_field
237  INTEGER(KIND=4), DIMENSION(mesh%me) :: i4_xml_field
238  INTEGER(KIND=1), DIMENSION(mesh%me) :: i1_xml_field
239  CHARACTER(LEN=200) :: ascii_or_binary
240 
241  IF (PRESENT(opt_st)) THEN
242  field_name=opt_st
243  ELSE
244  field_name='field'
245  END IF
246 
247  IF (SIZE(field)==0) RETURN
248 
249  IF (inputs%if_xml) THEN
250  ascii_or_binary = 'binary'
251  OPEN (unit=unit_file, file=file_name//'.vtu', status = 'unknown')
252  WRITE(unit_file,'(A)') '<?xml version="1.0" ?>'
253  WRITE(unit_file,'(A)', advance="no") '<VTKFile type="UnstructuredGrid" version="0.1" '
254  WRITE(unit_file,'(A)') 'compressor="vtkZLibDataCompressor" byte_order="LittleEndian">'
255  ELSE
256  ascii_or_binary = 'ascii'
257  OPEN (unit=unit_file, file=file_name//'.vtu',&
258  form = 'formatted', status = 'unknown')
259  WRITE(unit_file,'(A)') '<VTKFile type="UnstructuredGrid" version="0.1"'// &
260  ' byte_order="LittleEndian">'
261  END IF
262 
263  WRITE(unit_file,'(A)') '<UnstructuredGrid>'
264  WRITE(unit_file,'(A,I9,A,I9,A)') '<Piece NumberOfPoints="', mesh%np, &
265  '" NumberOfCells="', mesh%me, '">'
266 
267  !===PointData Block
268  WRITE(unit_file,'(A)') '<PointData>'
269  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="'//trim(adjustl(field_name))//&
270  '_cos" format="'//trim(adjustl(ascii_or_binary))//'" NumberOfComponents="3">'
271  IF (inputs%if_xml) THEN
272  DO n = 1, mesh%np
273  r4_threed_xml_field(3*(n-1)+1:3*n) = REAL(field(n,1:5:2),4)
274  END DO
275  CALL write_compressed_block(unit_file, r4_threed_xml_field(1:3*mesh%np))
276  ELSE
277  DO n = 1, mesh%np
278  WRITE(unit_file,'(3(e14.7,x))') field(n,1), field(n,2), field(n,3)
279  END DO
280  END IF
281  WRITE(unit_file,'(A)') '</DataArray>'
282  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="'//trim(adjustl(field_name))//&
283  '_sin" format="'//trim(adjustl(ascii_or_binary))//'" NumberOfComponents="3">'
284  IF (inputs%if_xml) THEN
285  DO n = 1, mesh%np
286  r4_threed_xml_field(3*(n-1)+1:3*n) = REAL(field(n,2:6:2),4)
287  END DO
288  CALL write_compressed_block(unit_file, r4_threed_xml_field(1:3*mesh%np))
289  ELSE
290  DO n = 1, mesh%np
291  WRITE(unit_file,'(3(e14.7,x))') field(n,2), field(n,4), field(n,6)
292  END DO
293  END IF
294  WRITE(unit_file,'(A)') '</DataArray>'
295  WRITE(unit_file,'(A)') '</PointData>'
296  !===End of PointData Block
297 
298  !===CellData Block
299  WRITE(unit_file,'(A)') '<CellData>'
300  WRITE(unit_file,'(A)') '</CellData>'
301  !===End of CellData Block
302 
303  !===Points Block
304  WRITE(unit_file,'(A)') '<Points>'
305  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="Points" '//&
306  'NumberOfComponents="3" format="'//trim(adjustl(ascii_or_binary))//'">'
307  IF (inputs%if_xml) THEN
308  DO n = 1, mesh%np
309  r4_threed_xml_field(3*(n-1)+1) = REAL(mesh%rr(1,n),4)
310  r4_threed_xml_field(3*(n-1)+2) = REAL(0.,4)
311  r4_threed_xml_field(3*(n-1)+3) = REAL(mesh%rr(2,n),4)
312  END DO
313  CALL write_compressed_block(unit_file, r4_threed_xml_field(1:3*mesh%np))
314  ELSE
315  DO n = 1, mesh%np
316  WRITE(unit_file,'(3(e14.7,x))') mesh%rr(1,n), 0.d0 , mesh%rr(3,n)
317  END DO
318  END IF
319  WRITE(unit_file,'(A)') '</DataArray>'
320  WRITE(unit_file,'(A)') '</Points>'
321  !===End of Points Block
322 
323  !===Cells Block
324  IF (mesh%gauss%n_w==3) THEN
325  type_cell = 5
326  ELSE IF (mesh%gauss%n_w==6) THEN
327  type_cell = 22
328  END IF
329  WRITE(unit_file,'(A)') '<Cells>'
330  WRITE(unit_file,'(A)') '<DataArray type="Int32" Name="connectivity" format="'&
331  //trim(adjustl(ascii_or_binary))//'">'
332  IF (inputs%if_xml) THEN
333  DO m = 1, mesh%me
334  i4_threed_xml_field(mesh%gauss%n_w*(m-1)+1) = int(mesh%jj(1,m)-1,4)
335  i4_threed_xml_field(mesh%gauss%n_w*(m-1)+2) = int(mesh%jj(2,m)-1,4)
336  i4_threed_xml_field(mesh%gauss%n_w*(m-1)+3) = int(mesh%jj(3,m)-1,4)
337  IF (type_cell==22) THEN
338  i4_threed_xml_field(mesh%gauss%n_w*(m-1)+4) = int(mesh%jj(6,m)-1,4)
339  i4_threed_xml_field(mesh%gauss%n_w*(m-1)+5) = int(mesh%jj(4,m)-1,4)
340  i4_threed_xml_field(mesh%gauss%n_w*(m-1)+6) = int(mesh%jj(5,m)-1,4)
341  END IF
342  END DO
343  CALL write_compressed_block(unit_file, i4_threed_xml_field)
344  ELSE
345  DO m = 1, mesh%me
346  WRITE(unit_file,'(3(I8,1x))') mesh%jj(1:3,m)-1
347  IF (type_cell==22) THEN
348  WRITE(unit_file,'(3(I8,1x))') mesh%jj(6,m)-1 , mesh%jj(4,m)-1 , mesh%jj(5,m)-1
349  END IF
350  END DO
351  END IF
352 
353  WRITE(unit_file,'(A)') '</DataArray>'
354  WRITE(unit_file,'(A)') '<DataArray type="Int32" Name="offsets" format="'&
355  //trim(adjustl(ascii_or_binary))//'">'
356  IF (inputs%if_xml) THEN
357  DO m = 1, mesh%me
358  i4_xml_field(m) = int(mesh%gauss%n_w*m,4)
359  END DO
360  CALL write_compressed_block(unit_file, i4_xml_field)
361  ELSE
362  DO m = 1, mesh%me
363  WRITE(unit_file,'(I8)') m*mesh%gauss%n_w
364  END DO
365  END IF
366 
367  WRITE(unit_file,'(A)') '</DataArray>'
368  WRITE(unit_file,'(A)') '<DataArray type="UInt8" Name="types" format="'&
369  //trim(adjustl(ascii_or_binary))//'">'
370  IF (inputs%if_xml) THEN
371  DO m = 1, mesh%me
372  i1_xml_field(m) = int(type_cell,1)
373  END DO
374  CALL write_compressed_block(unit_file, i1_xml_field)
375  ELSE
376  DO m = 1, mesh%me
377  WRITE(unit_file,'(I8)') type_cell
378  END DO
379  END IF
380  WRITE(unit_file,'(A)') '</DataArray>'
381  WRITE(unit_file,'(A)') '</Cells>'
382  ! End of Cells Block ---------------------------------------------------------
383 
384  WRITE(unit_file,'(A)') '</Piece>'
385  WRITE(unit_file,'(A)') '</UnstructuredGrid>'
386  WRITE(unit_file,'(A)') '</VTKFile>'
387 
388  CLOSE(unit_file)
389  END SUBROUTINE create_xml_vtu_vect_file
390 
391 
392 
393  SUBROUTINE make_vtu_file_arpack(communicator, mesh, header, field, field_name, what, num_vp)
394  USE def_type_mesh
395  USE my_util
396  IMPLICIT NONE
397  TYPE(mesh_type), INTENT(IN) :: mesh
398  CHARACTER(*), INTENT(IN) :: header
399  CHARACTER(*), INTENT(IN) :: field_name, what
400  INTEGER, INTENT(IN) :: num_vp
401  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: field
402  CHARACTER(LEN=200), DIMENSION(1) :: file_list
403  CHARACTER(LEN=3) :: st_rank
404 !#include "petsc/finclude/petsc.h"
405  petscerrorcode :: ierr
406  petscmpiint :: rank, nb_procs
407  mpi_comm :: communicator
408  CALL mpi_comm_rank(communicator, rank, ierr)
409  CALL mpi_comm_size(communicator, nb_procs, ierr)
410 
411  WRITE(st_rank,'(I3)') num_vp
412  file_list(1) = trim(header)//'_eigen_'//trim(adjustl(st_rank))
413 
414  CALL create_pvd_file(file_list, trim(header), num_vp, trim(what))
415 
416 !=========TEST FL Feb. 11th, 2013
417  IF (SIZE(field,2) == 6) THEN
418  !CALL create_vtu_vect_file(field, mesh, TRIM(ADJUSTL(file_list(1))), field_name)
419  CALL create_xml_vtu_vect_file(field, mesh, trim(adjustl(file_list(1))), field_name)
420  ELSE IF (SIZE(field,2) == 1) THEN
421  CALL create_xml_vtu_scal_file(field(:,1), mesh, trim(adjustl(file_list(1))), field_name)
422  ELSE IF (SIZE(field,2) .GT. 0) THEN
423  CALL create_xml_vtu_scal_file(field(:,1), mesh, trim(adjustl(file_list(1))), field_name)
424  ELSE
425  CALL error_petsc('Bug in make_vtu_file_arpack: field needs at least one component')
426  END IF
427 !=========TEST FL Feb. 11th, 2013
428  END SUBROUTINE make_vtu_file_arpack
429 
430  SUBROUTINE create_vtu_file_axi3d(field, mesh, file_name, opt_st)
431  USE def_type_mesh
432  IMPLICIT NONE
433  TYPE(mesh_type) :: mesh
434  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: field
435  CHARACTER(*), INTENT(IN) :: file_name
436  CHARACTER(*), OPTIONAL, INTENT(IN) :: opt_st
437  CHARACTER(LEN=200) :: field_name
438  REAL(KIND=8) :: theta, pi, dtheta
439  INTEGER :: unit_file=789, m, i, type_cell, &
440  nb_angle, k
441 
442  IF (SIZE(field,2)==0) RETURN
443 
444  IF (PRESENT(opt_st)) THEN
445  field_name=opt_st
446  ELSE
447  field_name='field'
448  END IF
449 
450  nb_angle = SIZE(field,1)
451  pi = acos(-1.d0)
452  dtheta = 2*pi/nb_angle
453 
454  OPEN (unit=unit_file, file=file_name//'.vtu',&
455  form = 'formatted', status = 'unknown')
456 
457  WRITE(unit_file,'(A)') '<VTKFile type="UnstructuredGrid" version="0.1"'// &
458  ' byte_order="LittleEndian">'
459  WRITE(unit_file,'(A)') '<UnstructuredGrid>'
460 
461  WRITE(unit_file,'(A,I9,A,I9,A)') '<Piece NumberOfPoints="', nb_angle*mesh%np, &
462  '" NumberOfCells="', nb_angle*mesh%me, '">'
463  ! PointData Block ------------------------------------------------------------
464  WRITE(unit_file,'(A)') '<PointData Scalars="truc">'
465 
466  IF (SIZE(field,2)==1) THEN
467  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="'//trim(adjustl(field_name))//&
468  '" format="ascii">'
469  DO k = 1, nb_angle
470  DO i = 1, mesh%np
471  WRITE(unit_file,'(e14.7)') field(k, 1, i)
472  ENDDO
473  END DO
474  ELSE
475  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="'//trim(adjustl(field_name))//&
476  '" format="ascii" NumberOfComponents="3">'
477  DO k = 1, nb_angle
478  DO i = 1, mesh%np
479  WRITE(unit_file,'(3(e14.7,x))') field(k, 1, i), field(k, 2, i), field(k, 3, i)
480  ENDDO
481  END DO
482  END IF
483 
484  WRITE(unit_file,'(A)') '</DataArray>'
485  WRITE(unit_file,'(A)') '</PointData>'
486  ! End of PointData Block -----------------------------------------------------
487 
488  ! CellData Block -------------------------------------------------------------
489  WRITE(unit_file,'(A)') '<CellData>'
490  WRITE(unit_file,'(A)') '</CellData>'
491  ! End of CellData Block ------------------------------------------------------
492 
493  ! Points Block ---------------------------------------------------------------
494  WRITE(unit_file,'(A)') '<Points>'
495  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="Points" '//&
496  'NumberOfComponents="3" format="ascii">'
497  DO k = 1, nb_angle
498  theta = (k-1)*dtheta
499  DO i=1, mesh%np
500  WRITE(unit_file,'(3(e14.7,2x))') mesh%rr(1,i)*cos(theta), &
501  mesh%rr(1,i)*sin(theta), mesh%rr(2,i)
502  ENDDO
503  END DO
504  WRITE(unit_file,'(A)') '</DataArray>'
505  WRITE(unit_file,'(A)') '</Points>'
506  ! End of Points Block --------------------------------------------------------
507 
508  ! Cells Block ----------------------------------------------------------------
509  type_cell = 13
510  WRITE(unit_file,'(A)') '<Cells>'
511  WRITE(unit_file,'(A)') '<DataArray type="Int64" Name="connectivity" format="ascii">'
512  DO k = 1, nb_angle-1
513  DO m = 1, mesh%me
514  WRITE(unit_file,'(3(I8,1x))') mesh%jj(1:3,m)-1+(k-1)*mesh%np
515  WRITE(unit_file,'(3(I8,1x))') mesh%jj(1:3,m)-1+k*mesh%np
516  END DO
517  END DO
518  k = nb_angle
519  DO m = 1, mesh%me
520  WRITE(unit_file,'(3(I8,1x))') mesh%jj(1:3,m)-1+(k-1)*mesh%np
521  WRITE(unit_file,'(3(I8,1x))') mesh%jj(1:3,m)-1
522  END DO
523  WRITE(unit_file,'(A)') '</DataArray>'
524  WRITE(unit_file,'(A)') '<DataArray type="Int64" Name="offsets" format="ascii">'
525  DO m = 1, nb_angle*mesh%me
526  WRITE(unit_file,'(I8)') 6*m
527  END DO
528  WRITE(unit_file,'(A)') '</DataArray>'
529  WRITE(unit_file,'(A)') '<DataArray type="UInt8" Name="types" format="ascii">'
530  DO m = 1, nb_angle*mesh%me
531  WRITE(unit_file,'(I8)') type_cell
532  END DO
533  WRITE(unit_file,'(A)') '</DataArray>'
534  WRITE(unit_file,'(A)') '</Cells>'
535  ! End of Cells Block ---------------------------------------------------------
536 
537  WRITE(unit_file,'(A)') '</Piece>'
538  WRITE(unit_file,'(A)') '</UnstructuredGrid>'
539  WRITE(unit_file,'(A)') '</VTKFile>'
540 
541  CLOSE(unit_file)
542  END SUBROUTINE create_vtu_file_axi3d
543 
544  SUBROUTINE make_vtu_file_3d(communicator, mesh, header, &
545  field, field_name, what, opt_it)
546  USE def_type_mesh
547  USE my_util
548  IMPLICIT NONE
549  TYPE(mesh_type), INTENT(IN) :: mesh
550  CHARACTER(*), INTENT(IN) :: header
551  CHARACTER(*), INTENT(IN) :: field_name, what
552  INTEGER, OPTIONAL, INTENT(IN) :: opt_it
553  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: field
554  INTEGER :: j, it
555  CHARACTER(LEN=200), DIMENSION(:), POINTER :: file_list
556  CHARACTER(LEN=3) :: st_rank, st_it
557 !#include "petsc/finclude/petsc.h"
558  petscerrorcode :: ierr
559  petscmpiint :: rank, nb_procs
560  mpi_comm :: communicator
561  CALL mpi_comm_rank(communicator, rank, ierr)
562  CALL mpi_comm_size(communicator, nb_procs, ierr)
563  ALLOCATE(file_list(nb_procs))
564 
565  IF (PRESENT(opt_it)) THEN
566  it = opt_it
567  WRITE(st_it,'(I3)') it
568  DO j = 1, nb_procs
569  WRITE(st_rank,'(I3)') j
570  file_list(j) = trim(header)//'_proc_'//trim(adjustl(st_rank))//&
571  '_it_'//trim(adjustl(st_it))
572  END DO
573  ELSE
574  DO j = 1, nb_procs
575  WRITE(st_rank,'(I3)') j
576  file_list(j) = trim(header)//'_proc_'//trim(adjustl(st_rank))
577  END DO
578  END IF
579 
580  CALL check_list(communicator, file_list, mesh%np)
581  IF (rank==0) THEN
582  IF (PRESENT(opt_it)) THEN
583  it = opt_it
584  ELSE
585  it = 1
586 
587  END IF
588  CALL create_pvd_file(file_list, trim(header), it, trim(what))
589  END IF
590 
591  !CALL create_vtu_file_3D(field, mesh, TRIM(ADJUSTL(file_list(rank+1))), &
592  ! opt_st=field_name)
593  CALL create_xml_vtu_file_3d(field, mesh, trim(adjustl(file_list(rank+1))), &
594  opt_st=field_name)
595 
596  END SUBROUTINE make_vtu_file_3d
597 
598 
599 
600  SUBROUTINE create_xml_vtu_file_3d(field, mesh, file_name, opt_st)
601  USE def_type_mesh
602  USE my_util
603  USE input_data
604  USE zlib_base64
605  IMPLICIT NONE
606  TYPE(mesh_type) :: mesh
607  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: field
608  CHARACTER(*), INTENT(IN) :: file_name
609  CHARACTER(*), OPTIONAL, INTENT(IN) :: opt_st
610  CHARACTER(LEN=200) :: field_name
611  INTEGER :: unit_file=789, m, n, type_cell, stride
612  REAL(KIND=4), DIMENSION(3*mesh%np) :: r4_threed_xml_field
613  INTEGER(KIND=4), DIMENSION(mesh%gauss%n_w*mesh%me) :: i4_threed_xml_field
614  INTEGER(KIND=4), DIMENSION(mesh%me) :: i4_xml_field
615  INTEGER(KIND=1), DIMENSION(mesh%me) :: i1_xml_field
616  CHARACTER(LEN=200) :: ascii_or_binary
617 
618  IF (SIZE(field,2)==0) RETURN
619 
620  IF (SIZE(mesh%jj,1)==6) THEN
621  type_cell = 13
622  stride = 6
623  ELSE IF (SIZE(mesh%jj,1)==15) THEN
624  type_cell = 26
625  stride = 15
626  ELSE
627  type_cell = 0
628  stride = 0
629  CALL error_petsc('Bug in create_vtu_file_3D: SIZE(mesh%jj,1) is wrong')
630  END IF
631 
632  IF (PRESENT(opt_st)) THEN
633  field_name=opt_st
634  ELSE
635  field_name='field'
636  END IF
637 
638  IF (inputs%if_xml) THEN
639  ascii_or_binary = 'binary'
640  OPEN (unit=unit_file, file=file_name//'.vtu', status = 'unknown')
641  WRITE(unit_file,'(A)') '<?xml version="1.0" ?>'
642  WRITE(unit_file,'(A)', advance="no") '<VTKFile type="UnstructuredGrid" version="0.1" '
643  WRITE(unit_file,'(A)') 'compressor="vtkZLibDataCompressor" byte_order="LittleEndian">'
644  ELSE
645  ascii_or_binary = 'ascii'
646  OPEN (unit=unit_file, file=file_name//'.vtu',&
647  form = 'formatted', status = 'unknown')
648  WRITE(unit_file,'(A)') '<VTKFile type="UnstructuredGrid" version="0.1"'// &
649  ' byte_order="LittleEndian">'
650  END IF
651 
652  WRITE(unit_file,'(A)') '<UnstructuredGrid>'
653  WRITE(unit_file,'(A,I9,A,I9,A)') '<Piece NumberOfPoints="', mesh%np, &
654  '" NumberOfCells="', mesh%me, '">'
655 
656  !===PointData Block
657  WRITE(unit_file,'(A)') '<PointData Scalars="truc">'
658  IF (SIZE(field,1)==1) THEN
659  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="'//trim(adjustl(field_name))//&
660  '" format="'//trim(adjustl(ascii_or_binary))//'">'
661  IF (inputs%if_xml) THEN
662  CALL write_compressed_block(unit_file, REAL(field(1,1:mesh%np),4))
663  ELSE
664  DO n = 1, mesh%np
665  WRITE(unit_file,'(e14.7)') field(1,n)
666  END DO
667  END IF
668  ELSE
669  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="'//trim(adjustl(field_name))//&
670  '" format="'//trim(adjustl(ascii_or_binary))//'" NumberOfComponents="3">'
671  IF (inputs%if_xml) THEN
672  DO n = 1, mesh%np
673  r4_threed_xml_field(3*(n-1)+1:3*n) = REAL(field(1:3,n),4)
674  END DO
675  CALL write_compressed_block(unit_file, r4_threed_xml_field(1:3*mesh%np))
676  ELSE
677  DO n = 1, mesh%np
678  WRITE(unit_file,'(3(e14.7,x))') field(1,n), field(2,n), field(3,n)
679  END DO
680  END IF
681  END IF
682  WRITE(unit_file,'(A)') '</DataArray>'
683  WRITE(unit_file,'(A)') '</PointData>'
684  !===End of PointData Block
685 
686  !===CellData Block
687  WRITE(unit_file,'(A)') '<CellData>'
688  WRITE(unit_file,'(A)') '</CellData>'
689  !===End of CellData Block
690 
691  !===Points Block
692  WRITE(unit_file,'(A)') '<Points>'
693  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="Points" '//&
694  'NumberOfComponents="3" format="'//trim(adjustl(ascii_or_binary))//'">'
695  IF (inputs%if_xml) THEN
696  DO n = 1, mesh%np
697  r4_threed_xml_field(3*(n-1)+1:3*n) = REAL(mesh%rr(1:3,n),4)
698  END DO
699  CALL write_compressed_block(unit_file, r4_threed_xml_field(1:3*mesh%np))
700  ELSE
701  DO n = 1, mesh%np
702  WRITE(unit_file,'(3(e14.7,x))') mesh%rr(1,n), mesh%rr(2,n), mesh%rr(3,n)
703  END DO
704  END IF
705  WRITE(unit_file,'(A)') '</DataArray>'
706  WRITE(unit_file,'(A)') '</Points>'
707  !===End of Points Block
708 
709  !===Cells Block
710  WRITE(unit_file,'(A)') '<Cells>'
711  WRITE(unit_file,'(A)') '<DataArray type="Int32" Name="connectivity" format="'&
712  //trim(adjustl(ascii_or_binary))//'">'
713  IF (inputs%if_xml) THEN
714  DO m = 1, mesh%me
715  i4_threed_xml_field(mesh%gauss%n_w*(m-1)+1:mesh%gauss%n_w*m) = int(mesh%jj(1:mesh%gauss%n_w,m)-1,4)
716  END DO
717  CALL write_compressed_block(unit_file, i4_threed_xml_field)
718  ELSE
719  DO m = 1, mesh%me
720  WRITE(unit_file,'(15(I8,1x))') mesh%jj(:,m)-1
721  END DO
722  END IF
723 
724  WRITE(unit_file,'(A)') '</DataArray>'
725  WRITE(unit_file,'(A)') '<DataArray type="Int32" Name="offsets" format="'&
726  //trim(adjustl(ascii_or_binary))//'">'
727  IF (inputs%if_xml) THEN
728  DO m = 1, mesh%me
729  i4_xml_field(m) = int(stride*m,4)
730  END DO
731  CALL write_compressed_block(unit_file, i4_xml_field)
732  ELSE
733  DO m = 1, mesh%me
734  WRITE(unit_file,'(I8)') stride*m
735  END DO
736  END IF
737 
738  WRITE(unit_file,'(A)') '</DataArray>'
739  WRITE(unit_file,'(A)') '<DataArray type="UInt8" Name="types" format="'&
740  //trim(adjustl(ascii_or_binary))//'">'
741  IF (inputs%if_xml) THEN
742  DO m = 1, mesh%me
743  i1_xml_field(m) = int(type_cell,1)
744  END DO
745  CALL write_compressed_block(unit_file, i1_xml_field)
746  ELSE
747  DO m = 1, mesh%me
748  WRITE(unit_file,'(I8)') type_cell
749  END DO
750  END IF
751  WRITE(unit_file,'(A)') '</DataArray>'
752  WRITE(unit_file,'(A)') '</Cells>'
753  !===End of Cells Block
754 
755  WRITE(unit_file,'(A)') '</Piece>'
756  WRITE(unit_file,'(A)') '</UnstructuredGrid>'
757  WRITE(unit_file,'(A)') '</VTKFile>'
758 
759  CLOSE(unit_file)
760  END SUBROUTINE create_xml_vtu_file_3d
761 
762 !!$ SUBROUTINE make_vtu_file_scalar_2D(comm, mesh, header, field, field_name, what, opt_it)
763 !!$ USE def_type_mesh
764 !!$ USE my_util
765 !!$ IMPLICIT NONE
766 !!$ TYPE(mesh_type), INTENT(IN) :: mesh
767 !!$ CHARACTER(*), INTENT(IN) :: header
768 !!$ CHARACTER(*), INTENT(IN) :: field_name, what
769 !!$ INTEGER, OPTIONAL, INTENT(IN) :: opt_it
770 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: field
771 !!$ INTEGER :: j, it
772 !!$ CHARACTER(LEN=200), DIMENSION(:), POINTER :: file_list
773 !!$ CHARACTER(LEN=3) :: st_rank, st_it
774 !!$!#include "petsc/finclude/petsc.h"
775 !!$ PetscErrorCode :: ierr
776 !!$ PetscMPIInt :: rank, nb_procs
777 !!$ MPI_Comm :: comm
778 !!$ CALL MPI_Comm_rank(comm, rank, ierr)
779 !!$ CALL MPI_Comm_size(comm, nb_procs, ierr)
780 !!$ ALLOCATE(file_list(nb_procs))
781 !!$ IF (PRESENT(opt_it)) THEN
782 !!$ it = opt_it
783 !!$ WRITE(st_it,'(I3)') it
784 !!$ DO j = 1, nb_procs
785 !!$ WRITE(st_rank,'(I3)') j
786 !!$ file_list(j) = TRIM(header)//'_proc_'//TRIM(ADJUSTL(st_rank))//&
787 !!$ '_it_'//TRIM(ADJUSTL(st_it))
788 !!$ END DO
789 !!$ ELSE
790 !!$ DO j = 1, nb_procs
791 !!$ WRITE(st_rank,'(I3)') j
792 !!$ file_list(j) = TRIM(header)//'_proc_'//TRIM(ADJUSTL(st_rank))
793 !!$ END DO
794 !!$ END IF
795 !!$
796 !!$ CALL check_list(comm, file_list, mesh%np)
797 !!$ IF (rank==0) THEN
798 !!$ IF (PRESENT(opt_it)) THEN
799 !!$ it = opt_it
800 !!$ ELSE
801 !!$ it = 1
802 !!$ END IF
803 !!$ CALL create_pvd_file(file_list, TRIM(header), it, TRIM(what))
804 !!$ END IF
805 !!$ CALL create_xml_vtu_scal_file(field, mesh, TRIM(ADJUSTL(file_list(rank+1))), TRIM(ADJUSTL(field_name)))
806 !!$ END SUBROUTINE make_vtu_file_scalar_2D
807 
808  !!$ SUBROUTINE create_vtk_file(comm, field, mesh, file_name, opt_it)
809 !!$ USE def_type_mesh
810 !!$ USE chaine_caractere
811 !!$ IMPLICIT NONE
812 !!$ TYPE(mesh_type) :: mesh
813 !!$ REAL(KIND=8), DIMENSION(:), INTENT(IN) :: field
814 !!$ CHARACTER(*), INTENT(IN) :: file_name
815 !!$ INTEGER, OPTIONAL :: opt_it
816 !!$ CHARACTER(len=4) :: tit
817 !!$ CHARACTER(len=3) :: st_it
818 !!$ INTEGER :: unit_file, l, lblank, start, i, j, nit
819 !!$ PetscErrorCode :: ierr
820 !!$ PetscMPIInt :: rank, nb_procs
821 !!$ MPI_Comm :: comm
822 !!$
823 !!$ CALL MPI_Comm_rank(comm,rank,ierr)
824 !!$ CALL MPI_Comm_size(comm,nb_procs,ierr)
825 !!$
826 !!$ IF (PRESENT(opt_it)) THEN
827 !!$ WRITE(st_it,'(i3)') opt_it
828 !!$ lblank = eval_blank(3,st_it)
829 !!$ DO l = 1, lblank - 1
830 !!$ st_it(l:l) = '0'
831 !!$ END DO
832 !!$ nit = opt_it
833 !!$ ELSE
834 !!$ st_it='001'
835 !!$ nit = 1
836 !!$ END IF
837 !!$
838 !!$ unit_file=rank+10
839 !!$ IF (rank==0) THEN
840 !!$ OPEN (UNIT=unit_file, FILE=file_name//'.visit', FORM = 'formatted', STATUS = 'unknown')
841 !!$ WRITE(unit_file,'(A,i7)') '!NBLOCKS ', nb_procs
842 !!$ DO j = 0, nb_procs-1
843 !!$ WRITE(tit,'(I4)') j
844 !!$ lblank = eval_blank(4,tit)
845 !!$ DO l = 1, lblank - 1
846 !!$ tit(l:l) = '0'
847 !!$ END DO
848 !!$ start = 4 - nb_digit() + 1
849 !!$ WRITE(unit_file,'(A)') file_name//'_'//tit(start:)//'_'//st_it//'.vtk'
850 !!$ END DO
851 !!$ CLOSE(unit_file) !Ecriture pour visit
852 !!$
853 !!$ OPEN (UNIT=unit_file, FILE=file_name//'.pvd', FORM = 'formatted', STATUS = 'unknown')
854 !!$ IF (nit == 1) THEN
855 !!$ WRITE(unit_file, '(A)') '<?xml version="1.0"?>'
856 !!$ WRITE(unit_file, '(A)') '<VTKFile type="Collection" version="0.1" '// &
857 !!$ 'byte_order="BigEndian" compressor="vtkZLibDataCompressor">'
858 !!$ WRITE(unit_file, '(A)') '<Collection>'
859 !!$ ELSE
860 !!$ DO j = 1, 3+(nit-1)*nb_procs
861 !!$ READ(unit_file,*)
862 !!$ END DO
863 !!$ END IF
864 !!$ DO j = 0, nb_procs-1
865 !!$ WRITE(tit,'(I4)') j
866 !!$ lblank = eval_blank(4,tit)
867 !!$ DO l = 1, lblank - 1
868 !!$ tit(l:l) = '0'
869 !!$ END DO
870 !!$ start = 4 - nb_digit() + 1
871 !!$ WRITE(unit_file,'(A)') '<DataSet timestep="'//st_it//'" group="" '// &
872 !!$ 'part="'//tit(2:)//'" file="./'//file_name//'_'//tit(start:)//'_'//st_it//'.vtu'//'"/>'
873 !!$ END DO
874 !!$ WRITE(unit_file, '(A)') '</Collection>'
875 !!$ WRITE(unit_file, '(A)') '</VTKFile>'
876 !!$
877 !!$ CLOSE(unit_file) !Ecriture pour paraview
878 !!$ END IF
879 !!$
880 !!$ IF (SIZE(field,1)==0) RETURN
881 !!$
882 !!$ WRITE(tit,'(I4)') rank
883 !!$ lblank = eval_blank(4,tit)
884 !!$ DO l = 1, lblank - 1
885 !!$ tit(l:l) = '0'
886 !!$ END DO
887 !!$ start = 4 - nb_digit() + 1
888 !!$
889 !!$ unit_file=rank+10
890 !!$ OPEN (UNIT=unit_file, FILE=file_name//'_'//tit(start:)//'_'//st_it//'.vtk', &
891 !!$ FORM = 'formatted', STATUS = 'unknown')
892 !!$ WRITE(unit_file,'(A)') '# vtk DataFile Version 3.0'
893 !!$ WRITE(unit_file,'(A)') 'vtk '//file_name//''
894 !!$ WRITE(unit_file,'(A)')'ASCII'
895 !!$ WRITE(unit_file,'(A)')'DATASET UNSTRUCTURED_GRID'
896 !!$
897 !!$ WRITE(unit_file,'(A,I7,A)')'POINTS ', mesh%np, ' float'
898 !!$ WRITE(*,*) 'points ...'
899 !!$ DO i=1, mesh%np
900 !!$ WRITE(unit_file,'(2(e14.7,2x),A)') mesh%rr(1,i), &
901 !!$ mesh%rr(2,i), ' 0.0 '
902 !!$ ENDDO
903 !!$ WRITE(*,*) 'cells ...'
904 !!$ IF (mesh%gauss%n_w==3) THEN
905 !!$ WRITE(unit_file,'(A,I7,I8)') 'CELLS ', mesh%me, mesh%me*4
906 !!$ DO i=1, mesh%me
907 !!$ WRITE(unit_file,'(A,3(I8,1x))') '3 ', mesh%jj(1,i)-1, mesh%jj(2,i)-1, mesh%jj(3,i)-1
908 !!$ ENDDO
909 !!$ WRITE(unit_file,'(A,I7)') 'CELL_TYPES ', mesh%me
910 !!$ DO i=1, mesh%me
911 !!$ WRITE(unit_file,'(A)') '5'
912 !!$ ENDDO
913 !!$ ELSE IF (mesh%gauss%n_w==6) THEN
914 !!$ WRITE(unit_file,'(A,I7,I8)') 'CELLS ', 4*mesh%me, 4*mesh%me*4
915 !!$ DO i=1, mesh%me
916 !!$ WRITE(unit_file,'(A,3(I8,1x))') '3 ', mesh%jj(1,i)-1, mesh%jj(6,i)-1, mesh%jj(5,i)-1
917 !!$ WRITE(unit_file,'(A,3(I8,1x))') '3 ', mesh%jj(2,i)-1, mesh%jj(4,i)-1, mesh%jj(6,i)-1
918 !!$ WRITE(unit_file,'(A,3(I8,1x))') '3 ', mesh%jj(3,i)-1, mesh%jj(5,i)-1, mesh%jj(4,i)-1
919 !!$ WRITE(unit_file,'(A,3(I8,1x))') '3 ', mesh%jj(6,i)-1, mesh%jj(4,i)-1, mesh%jj(5,i)-1
920 !!$ ENDDO
921 !!$ WRITE(unit_file,'(A,I7)') 'CELL_TYPES ', 4*mesh%me
922 !!$
923 !!$ DO i=1, 4*mesh%me
924 !!$ WRITE(unit_file,'(A)') '5'
925 !!$ ENDDO
926 !!$ END IF
927 !!$
928 !!$ WRITE(*,*) 'data ...'
929 !!$ WRITE(unit_file,'(A,I7)') 'POINT_DATA ',mesh%np
930 !!$ WRITE(unit_file,'(A)') 'SCALARS scalars float 1'
931 !!$ WRITE(unit_file,'(A)') 'LOOKUP_TABLE default'
932 !!$ DO i=1, mesh%np
933 !!$ WRITE(unit_file,'(e14.7,2x)') field(i)
934 !!$ ENDDO
935 !!$ CLOSE(unit_file)
936 !!$ END SUBROUTINE create_vtk_file
937 
938 !!$ SUBROUTINE create_vtu_vect_file(field, mesh, file_name, opt_st)
939 !!$ USE def_type_mesh
940 !!$ IMPLICIT NONE
941 !!$ TYPE(mesh_type) :: mesh
942 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: field
943 !!$ CHARACTER(*), INTENT(IN) :: file_name
944 !!$ INTEGER :: unit_file=789, m, i, type_cell
945 !!$ CHARACTER(*), OPTIONAL, INTENT(IN) :: opt_st
946 !!$ CHARACTER(LEN=200) :: field_name
947 !!$
948 !!$ IF (PRESENT(opt_st)) THEN
949 !!$ field_name=opt_st
950 !!$ ELSE
951 !!$ field_name='field'
952 !!$ END IF
953 !!$
954 !!$ IF (SIZE(field)==0) RETURN
955 !!$ OPEN (UNIT=unit_file, FILE=file_name//'.vtu',&
956 !!$ FORM = 'formatted', STATUS = 'unknown')
957 !!$
958 !!$ WRITE(unit_file,'(A)') '<VTKFile type="UnstructuredGrid" version="0.1"'// &
959 !!$ ' byte_order="LittleEndian">'
960 !!$ WRITE(unit_file,'(A)') '<UnstructuredGrid>'
961 !!$
962 !!$ WRITE(unit_file,'(A,I9,A,I9,A)') '<Piece NumberOfPoints="', mesh%np, &
963 !!$ '" NumberOfCells="', mesh%me, '">'
964 !!$ ! PointData Block ------------------------------------------------------------
965 !!$ WRITE(unit_file,'(A)') '<PointData>'
966 !!$ WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="'//TRIM(ADJUSTL(field_name))//&
967 !!$ '_cos" format="ascii" NumberOfComponents="3">'
968 !!$ DO i = 1, mesh%np
969 !!$ WRITE(unit_file,'(e14.7)') field(i,1), field(i,3), field(i,5)
970 !!$ ENDDO
971 !!$ WRITE(unit_file,'(A)') '</DataArray>'
972 !!$ WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="'//TRIM(ADJUSTL(field_name))//&
973 !!$ '_sin" format="ascii" NumberOfComponents="3">'
974 !!$ DO i = 1, mesh%np
975 !!$ WRITE(unit_file,'(e14.7)') field(i,2), field(i,4), field(i,6)
976 !!$ ENDDO
977 !!$ WRITE(unit_file,'(A)') '</DataArray>'
978 !!$ WRITE(unit_file,'(A)') '</PointData>'
979 !!$ ! End of PointData Block -----------------------------------------------------
980 !!$
981 !!$ ! CellData Block -------------------------------------------------------------
982 !!$ WRITE(unit_file,'(A)') '<CellData>'
983 !!$ WRITE(unit_file,'(A)') '</CellData>'
984 !!$ ! End of CellData Block ------------------------------------------------------
985 !!$
986 !!$ ! Points Block ---------------------------------------------------------------
987 !!$ WRITE(unit_file,'(A)') '<Points>'
988 !!$ WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="Points" '//&
989 !!$ 'NumberOfComponents="3" format="ascii">'
990 !!$ DO i=1, mesh%np
991 !!$ WRITE(unit_file,'(e14.7,A,e14.7)') mesh%rr(1,i), ' 0.0 ' , &
992 !!$ mesh%rr(2,i)
993 !!$ ENDDO
994 !!$ WRITE(unit_file,'(A)') '</DataArray>'
995 !!$ WRITE(unit_file,'(A)') '</Points>'
996 !!$ ! End of Points Block --------------------------------------------------------
997 !!$
998 !!$ ! Cells Block ----------------------------------------------------------------
999 !!$ IF (mesh%gauss%n_w==3) THEN
1000 !!$ type_cell = 5
1001 !!$ ELSE IF (mesh%gauss%n_w==6) THEN
1002 !!$ type_cell = 22
1003 !!$ END IF
1004 !!$ WRITE(unit_file,'(A)') '<Cells>'
1005 !!$ WRITE(unit_file,'(A)') '<DataArray type="Int64" Name="connectivity" format="ascii">'
1006 !!$ DO m = 1, mesh%me
1007 !!$ WRITE(unit_file,'(3(I8,1x))') mesh%jj(1:3,m)-1
1008 !!$ IF (type_cell==22) THEN
1009 !!$ WRITE(unit_file,'(3(I8,1x))') mesh%jj(6,m)-1 , mesh%jj(4,m)-1 , mesh%jj(5,m)-1
1010 !!$ END IF
1011 !!$ END DO
1012 !!$ WRITE(unit_file,'(A)') '</DataArray>'
1013 !!$ WRITE(unit_file,'(A)') '<DataArray type="Int64" Name="offsets" format="ascii">'
1014 !!$ DO m = 1, mesh%me
1015 !!$ WRITE(unit_file,'(I8)') m*mesh%gauss%n_w
1016 !!$ END DO
1017 !!$ WRITE(unit_file,'(A)') '</DataArray>'
1018 !!$ WRITE(unit_file,'(A)') '<DataArray type="UInt8" Name="types" format="ascii">'
1019 !!$ DO m = 1, mesh%me
1020 !!$ WRITE(unit_file,'(I8)') type_cell
1021 !!$ END DO
1022 !!$ WRITE(unit_file,'(A)') '</DataArray>'
1023 !!$ WRITE(unit_file,'(A)') '</Cells>'
1024 !!$ ! End of Cells Block ---------------------------------------------------------
1025 !!$
1026 !!$ WRITE(unit_file,'(A)') '</Piece>'
1027 !!$ WRITE(unit_file,'(A)') '</UnstructuredGrid>'
1028 !!$ WRITE(unit_file,'(A)') '</VTKFile>'
1029 !!$
1030 !!$ CLOSE(unit_file)
1031 !!$ END SUBROUTINE create_vtu_vect_file
1032 
1033  !!$ SUBROUTINE make_vtu_file_axi3D(communicator, mesh, header, &
1034 !!$ field, field_name, what, opt_it)
1035 !!$ USE def_type_mesh
1036 !!$ USE my_util
1037 !!$ IMPLICIT NONE
1038 !!$ TYPE(mesh_type), INTENT(IN) :: mesh
1039 !!$ CHARACTER(*), INTENT(IN) :: header
1040 !!$ CHARACTER(*), INTENT(IN) :: field_name, what
1041 !!$ INTEGER, OPTIONAL, INTENT(IN) :: opt_it
1042 !!$ REAL(KIND=8), DIMENSION(:,:,:),INTENT(IN) :: field
1043 !!$ INTEGER :: j, it
1044 !!$ CHARACTER(LEN=200), DIMENSION(:), POINTER :: file_list
1045 !!$ CHARACTER(LEN=3) :: st_rank
1046 !!$!#include "petsc/finclude/petsc.h"
1047 !!$ PetscErrorCode :: ierr
1048 !!$ PetscMPIInt :: rank, nb_procs
1049 !!$ MPI_Comm :: communicator
1050 !!$ CALL MPI_Comm_rank(communicator, rank, ierr)
1051 !!$ CALL MPI_Comm_size(communicator, nb_procs, ierr)
1052 !!$ ALLOCATE(file_list(nb_procs))
1053 !!$ DO j = 1, nb_procs
1054 !!$ WRITE(st_rank,'(I3)') j
1055 !!$ file_list(j) = TRIM(header)//'_proc_'//TRIM(ADJUSTL(st_rank))
1056 !!$ END DO
1057 !!$ CALL check_list(communicator, file_list, mesh%np)
1058 !!$ IF (rank==0) THEN
1059 !!$ IF (PRESENT(opt_it)) THEN
1060 !!$ it = opt_it
1061 !!$ ELSE
1062 !!$ it = 1
1063 !!$
1064 !!$ END IF
1065 !!$ CALL create_pvd_file(file_list, TRIM(header), it, TRIM(what))
1066 !!$ END IF
1067 !!$
1068 !!$ CALL create_vtu_file_axi3D(field, mesh, TRIM(ADJUSTL(file_list(rank+1))), &
1069 !!$ opt_st=field_name)
1070 !!$ END SUBROUTINE make_vtu_file_axi3D
1071 
1072 !!$ SUBROUTINE create_vtu_file_3D(field, mesh, file_name, opt_st)
1073 !!$ USE def_type_mesh
1074 !!$ USE my_util
1075 !!$ IMPLICIT NONE
1076 !!$ TYPE(mesh_type) :: mesh
1077 !!$ REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: field
1078 !!$ CHARACTER(*), INTENT(IN) :: file_name
1079 !!$ CHARACTER(*), OPTIONAL, INTENT(IN) :: opt_st
1080 !!$ CHARACTER(LEN=200) :: field_name
1081 !!$ INTEGER :: unit_file=789, m, n, type_cell, stride
1082 !!$
1083 !!$ IF (SIZE(field,2)==0) RETURN
1084 !!$
1085 !!$ IF (SIZE(mesh%jj,1)==6) THEN
1086 !!$ type_cell = 13
1087 !!$ stride = 6
1088 !!$ ELSE IF (SIZE(mesh%jj,1)==15) THEN
1089 !!$ type_cell = 26
1090 !!$ stride = 15
1091 !!$ ELSE
1092 !!$ type_cell = 0
1093 !!$ stride = 0
1094 !!$ CALL error_petsc('Bug in create_vtu_file_3D: SIZE(mesh%jj,1) is wrong')
1095 !!$ END IF
1096 !!$
1097 !!$ IF (PRESENT(opt_st)) THEN
1098 !!$ field_name=opt_st
1099 !!$ ELSE
1100 !!$ field_name='field'
1101 !!$ END IF
1102 !!$
1103 !!$ OPEN (UNIT=unit_file, FILE=file_name//'.vtu',&
1104 !!$ FORM = 'formatted', STATUS = 'unknown')
1105 !!$
1106 !!$ WRITE(unit_file,'(A)') '<VTKFile type="UnstructuredGrid" version="0.1"'// &
1107 !!$ ' byte_order="LittleEndian">'
1108 !!$ WRITE(unit_file,'(A)') '<UnstructuredGrid>'
1109 !!$
1110 !!$ WRITE(unit_file,'(A,I9,A,I9,A)') '<Piece NumberOfPoints="', mesh%np, &
1111 !!$ '" NumberOfCells="', mesh%me, '">'
1112 !!$ ! PointData Block ------------------------------------------------------------
1113 !!$ WRITE(unit_file,'(A)') '<PointData Scalars="truc">'
1114 !!$
1115 !!$ IF (SIZE(field,1)==1) THEN
1116 !!$ WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="'//TRIM(ADJUSTL(field_name))//&
1117 !!$ '" format="ascii">'
1118 !!$ DO n = 1, mesh%np
1119 !!$ WRITE(unit_file,'(e14.7)') field(1,n)
1120 !!$ END DO
1121 !!$ ELSE
1122 !!$ WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="'//TRIM(ADJUSTL(field_name))//&
1123 !!$ '" format="ascii" NumberOfComponents="3">'
1124 !!$ DO n = 1, mesh%np
1125 !!$ WRITE(unit_file,'(3(e14.7,x))') field(1,n), field(2,n), field(3,n)
1126 !!$ END DO
1127 !!$ END IF
1128 !!$
1129 !!$ WRITE(unit_file,'(A)') '</DataArray>'
1130 !!$ WRITE(unit_file,'(A)') '</PointData>'
1131 !!$ ! End of PointData Block -----------------------------------------------------
1132 !!$
1133 !!$ ! CellData Block -------------------------------------------------------------
1134 !!$ WRITE(unit_file,'(A)') '<CellData>'
1135 !!$ WRITE(unit_file,'(A)') '</CellData>'
1136 !!$ ! End of CellData Block ------------------------------------------------------
1137 !!$
1138 !!$ ! Points Block ---------------------------------------------------------------
1139 !!$ WRITE(unit_file,'(A)') '<Points>'
1140 !!$ WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="Points" '//&
1141 !!$ 'NumberOfComponents="3" format="ascii">'
1142 !!$ DO n = 1, mesh%np
1143 !!$ WRITE(unit_file,'(3(e14.7,x))') mesh%rr(1,n), mesh%rr(2,n), mesh%rr(3,n)
1144 !!$ END DO
1145 !!$ WRITE(unit_file,'(A)') '</DataArray>'
1146 !!$ WRITE(unit_file,'(A)') '</Points>'
1147 !!$ ! End of Points Block --------------------------------------------------------
1148 !!$
1149 !!$ ! Cells Block ----------------------------------------------------------------
1150 !!$ WRITE(unit_file,'(A)') '<Cells>'
1151 !!$ WRITE(unit_file,'(A)') '<DataArray type="Int64" Name="connectivity" format="ascii">'
1152 !!$ DO m = 1, mesh%me
1153 !!$ WRITE(unit_file,'(15(I8,1x))') mesh%jj(:,m)-1
1154 !!$ END DO
1155 !!$ WRITE(unit_file,'(A)') '</DataArray>'
1156 !!$ WRITE(unit_file,'(A)') '<DataArray type="Int64" Name="offsets" format="ascii">'
1157 !!$ DO m = 1, mesh%me
1158 !!$ WRITE(unit_file,'(I8)') stride*m
1159 !!$ END DO
1160 !!$ WRITE(unit_file,'(A)') '</DataArray>'
1161 !!$ WRITE(unit_file,'(A)') '<DataArray type="UInt8" Name="types" format="ascii">'
1162 !!$ DO m = 1, mesh%me
1163 !!$ WRITE(unit_file,'(I8)') type_cell
1164 !!$ END DO
1165 !!$ WRITE(unit_file,'(A)') '</DataArray>'
1166 !!$ WRITE(unit_file,'(A)') '</Cells>'
1167 !!$ ! End of Cells Block ---------------------------------------------------------
1168 !!$
1169 !!$ WRITE(unit_file,'(A)') '</Piece>'
1170 !!$ WRITE(unit_file,'(A)') '</UnstructuredGrid>'
1171 !!$ WRITE(unit_file,'(A)') '</VTKFile>'
1172 !!$
1173 !!$ CLOSE(unit_file)
1174 !!$ END SUBROUTINE create_vtu_file_3D
1175 
1176 !!$ FUNCTION nb_digit() RESULT(dg)
1177 !!$ USE my_util
1178 !!$ IMPLICIT NONE
1179 !!$ INTEGER :: code, nb_procs, dg
1180 !!$!#include "petsc/finclude/petsc.h"
1181 !!$ CALL MPI_Comm_size(PETSC_COMM_WORLD,nb_procs,code)
1182 !!$
1183 !!$ IF (nb_procs>9999) THEN
1184 !!$ CALL error_Petsc('nb_procs>9999')
1185 !!$ END IF
1186 !!$
1187 !!$ IF (nb_procs < 10) THEN
1188 !!$ dg=1
1189 !!$ ELSE IF (nb_procs < 99) THEN
1190 !!$ dg=2
1191 !!$ ELSE IF (nb_procs < 999) THEN
1192 !!$ dg=3
1193 !!$ ELSE IF (nb_procs < 9999) THEN
1194 !!$ dg=4
1195 !!$ ELSE
1196 !!$ dg=0
1197 !!$ END IF
1198 !!$ END FUNCTION nb_digit
1199 !!$
1200 END MODULE vtk_viz
vtk_viz::create_vtu_file_axi3d
subroutine, public create_vtu_file_axi3d(field, mesh, file_name, opt_st)
Definition: plot_vtk.f90:431
vtk_viz::check_list
subroutine, public check_list(communicator, file_list, check)
Definition: plot_vtk.f90:14
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:74
my_util
Definition: my_util.f90:1
my_util::error_petsc
subroutine error_petsc(string)
Definition: my_util.f90:16
input_data::inputs
type(my_data), public inputs
Definition: read_my_data.f90:237
vtk_viz::make_vtu_file_arpack
subroutine, public make_vtu_file_arpack(communicator, mesh, header, field, field_name, what, num_vp)
Definition: plot_vtk.f90:394
vtk_viz
Definition: plot_vtk.f90:1
vtk_viz::create_xml_vtu_scal_file
subroutine, public create_xml_vtu_scal_file(field, mesh, file_name, field_name)
Definition: plot_vtk.f90:82
vtk_viz::create_xml_vtu_file_3d
subroutine create_xml_vtu_file_3d(field, mesh, file_name, opt_st)
Definition: plot_vtk.f90:601
vtk_viz::create_pvd_file
subroutine, public create_pvd_file(file_list, file_header, time_step, what)
Definition: plot_vtk.f90:47
input_data
Definition: read_my_data.f90:233
vtk_viz::make_vtu_file_3d
subroutine, public make_vtu_file_3d(communicator, mesh, header, field, field_name, what, opt_it)
Definition: plot_vtk.f90:546
def_type_mesh
Definition: def_type_mesh.f90:4
vtk_viz::create_xml_vtu_vect_file
subroutine, public create_xml_vtu_vect_file(field, mesh, file_name, opt_st)
Definition: plot_vtk.f90:225
form
section doc_intro_frame_work_num_app Numerical approximation subsection doc_intro_fram_work_num_app_Fourier_FEM Fourier Finite element representation The SFEMaNS code uses a hybrid Fourier Finite element formulation The Fourier decomposition allows to approximate the problem’s solutions for each Fourier mode modulo nonlinear terms that are made explicit The variables are then approximated on a meridian section of the domain with a finite element method The numerical approximation of a function f $f f is written in the following generic form
Definition: doc_intro.h:190