SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
plot_vtk.f90
Go to the documentation of this file.
1 MODULE vtk_viz
2 
3  PUBLIC :: create_vtk_file, create_vtu_scal_file, create_pvd_file, check_list,&
4  create_vtu_vect_file, make_vtu_file_2d, make_vtu_file_arpack
5  PUBLIC :: create_vtu_file_axi3d
6  PUBLIC :: make_vtu_file_axi3d
7  PUBLIC :: create_vtu_file_3d
8  PUBLIC :: make_vtu_file_3d
9  PUBLIC :: make_vtu_file_scalar_2d
10  PRIVATE
11 
12 CONTAINS
13 
14  SUBROUTINE create_vtk_file(comm, field, mesh, file_name, opt_it)
15  USE def_type_mesh
16  USE chaine_caractere
17  IMPLICIT NONE
18  TYPE(mesh_type) :: mesh
19  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: field
20  CHARACTER(*), INTENT(IN) :: file_name
21  INTEGER, OPTIONAL :: opt_it
22  CHARACTER(len=4) :: tit
23  CHARACTER(len=3) :: st_it
24  INTEGER :: unit_file, l, lblank, start, i, j, nit
25 #include "petsc/finclude/petsc.h"
26  petscerrorcode :: ierr
27  petscmpiint :: rank, nb_procs
28  mpi_comm :: comm
29 
30  CALL mpi_comm_rank(comm,rank,ierr)
31  CALL mpi_comm_size(comm,nb_procs,ierr)
32 
33  IF (present(opt_it)) THEN
34  WRITE(st_it,'(i3)') opt_it
35  lblank = eval_blank(3,st_it)
36  DO l = 1, lblank - 1
37  st_it(l:l) = '0'
38  END DO
39  nit = opt_it
40  ELSE
41  st_it='001'
42  nit = 1
43  END IF
44 
45  unit_file=rank+10
46  IF (rank==0) THEN
47  OPEN (unit=unit_file, file=file_name//'.visit', form = 'formatted', status = 'unknown')
48  WRITE(unit_file,'(A,i7)') '!NBLOCKS ', nb_procs
49  DO j = 0, nb_procs-1
50  WRITE(tit,'(I4)') j
51  lblank = eval_blank(4,tit)
52  DO l = 1, lblank - 1
53  tit(l:l) = '0'
54  END DO
55  start = 4 - nb_digit() + 1
56  WRITE(unit_file,'(A)') file_name//'_'//tit(start:)//'_'//st_it//'.vtk'
57  END DO
58  CLOSE(unit_file) !Ecriture pour visit
59 
60  OPEN (unit=unit_file, file=file_name//'.pvd', form = 'formatted', status = 'unknown')
61  IF (nit == 1) THEN
62  WRITE(unit_file, '(A)') '<?xml version="1.0"?>'
63  WRITE(unit_file, '(A)') '<VTKFile type="Collection" version="0.1" '// &
64  'byte_order="BigEndian" compressor="vtkZLibDataCompressor">'
65  WRITE(unit_file, '(A)') '<Collection>'
66  ELSE
67  DO j = 1, 3+(nit-1)*nb_procs
68  READ(unit_file,*)
69  END DO
70  END IF
71  DO j = 0, nb_procs-1
72  WRITE(tit,'(I4)') j
73  lblank = eval_blank(4,tit)
74  DO l = 1, lblank - 1
75  tit(l:l) = '0'
76  END DO
77  start = 4 - nb_digit() + 1
78  WRITE(unit_file,'(A)') '<DataSet timestep="'//st_it//'" group="" '// &
79  'part="'//tit(2:)//'" file="./'//file_name//'_'//tit(start:)//'_'//st_it//'.vtu'//'"/>'
80  END DO
81  WRITE(unit_file, '(A)') '</Collection>'
82  WRITE(unit_file, '(A)') '</VTKFile>'
83 
84  CLOSE(unit_file) !Ecriture pour paraview
85  END IF
86 
87  IF (SIZE(field,1)==0) RETURN
88 
89  WRITE(tit,'(I4)') rank
90  lblank = eval_blank(4,tit)
91  DO l = 1, lblank - 1
92  tit(l:l) = '0'
93  END DO
94  start = 4 - nb_digit() + 1
95 
96  unit_file=rank+10
97  OPEN (unit=unit_file, file=file_name//'_'//tit(start:)//'_'//st_it//'.vtk', &
98  form = 'formatted', status = 'unknown')
99  WRITE(unit_file,'(A)') '# vtk DataFile Version 3.0'
100  WRITE(unit_file,'(A)') 'vtk '//file_name//''
101  WRITE(unit_file,'(A)')'ASCII'
102  WRITE(unit_file,'(A)')'DATASET UNSTRUCTURED_GRID'
103 
104  WRITE(unit_file,'(A,I7,A)')'POINTS ', mesh%np, ' float'
105  WRITE(*,*) 'points ...'
106  DO i=1, mesh%np
107  WRITE(unit_file,'(2(e14.7,2x),A)') mesh%rr(1,i), &
108  mesh%rr(2,i), ' 0.0 '
109  ENDDO
110  WRITE(*,*) 'cells ...'
111  IF (mesh%gauss%n_w==3) THEN
112  WRITE(unit_file,'(A,I7,I8)') 'CELLS ', mesh%me, mesh%me*4
113  DO i=1, mesh%me
114  WRITE(unit_file,'(A,3(I8,1x))') '3 ', mesh%jj(1,i)-1, mesh%jj(2,i)-1, mesh%jj(3,i)-1
115  ENDDO
116  WRITE(unit_file,'(A,I7)') 'CELL_TYPES ', mesh%me
117  DO i=1, mesh%me
118  WRITE(unit_file,'(A)') '5'
119  ENDDO
120  ELSE IF (mesh%gauss%n_w==6) THEN
121  WRITE(unit_file,'(A,I7,I8)') 'CELLS ', 4*mesh%me, 4*mesh%me*4
122  DO i=1, mesh%me
123  WRITE(unit_file,'(A,3(I8,1x))') '3 ', mesh%jj(1,i)-1, mesh%jj(6,i)-1, mesh%jj(5,i)-1
124  WRITE(unit_file,'(A,3(I8,1x))') '3 ', mesh%jj(2,i)-1, mesh%jj(4,i)-1, mesh%jj(6,i)-1
125  WRITE(unit_file,'(A,3(I8,1x))') '3 ', mesh%jj(3,i)-1, mesh%jj(5,i)-1, mesh%jj(4,i)-1
126  WRITE(unit_file,'(A,3(I8,1x))') '3 ', mesh%jj(6,i)-1, mesh%jj(4,i)-1, mesh%jj(5,i)-1
127  ENDDO
128  WRITE(unit_file,'(A,I7)') 'CELL_TYPES ', 4*mesh%me
129 
130  DO i=1, 4*mesh%me
131  WRITE(unit_file,'(A)') '5'
132  ENDDO
133  END IF
134 
135  WRITE(*,*) 'data ...'
136  WRITE(unit_file,'(A,I7)') 'POINT_DATA ',mesh%np
137  WRITE(unit_file,'(A)') 'SCALARS scalars float 1'
138  WRITE(unit_file,'(A)') 'LOOKUP_TABLE default'
139  DO i=1, mesh%np
140  WRITE(unit_file,'(e14.7,2x)') field(i)
141  ENDDO
142  CLOSE(unit_file)
143  END SUBROUTINE create_vtk_file
144 
145  SUBROUTINE check_list(communicator, file_list, check)
146  IMPLICIT NONE
147  CHARACTER(LEN=200), DIMENSION(:), POINTER :: file_list
148  CHARACTER(LEN=200), DIMENSION(:), POINTER :: dummy_list
149  INTEGER, DIMENSION(SIZE(file_list)) :: check_mylist
150  INTEGER :: check, n, count
151 #include "petsc/finclude/petsc.h"
152  mpi_comm :: communicator
153  petscmpiint :: rank, nb_procs
154  petscerrorcode :: ierr
155  CALL mpi_comm_rank(communicator, rank, ierr)
156  CALL mpi_comm_size(communicator, nb_procs, ierr)
157 
158  CALL mpi_allgather(check, 1, mpi_integer, check_mylist, 1, &
159  mpi_integer, communicator, ierr)
160 
161  count = 0
162  DO n = 1, SIZE(file_list)
163  IF (check_mylist(n)==0) cycle
164  count = count + 1
165  END DO
166  ALLOCATE(dummy_list(count))
167  count = 0
168  DO n = 1, SIZE(file_list)
169  IF (check_mylist(n)==0) cycle
170  count = count + 1
171  dummy_list(count) = file_list(n)
172  END DO
173  DEALLOCATE(file_list)
174  ALLOCATE(file_list(count))
175  file_list = dummy_list
176  END SUBROUTINE check_list
177 
178  SUBROUTINE create_pvd_file(file_list, file_header, time_step, what)
179  IMPLICIT NONE
180  CHARACTER(*), DIMENSION(:), INTENT(IN) :: file_list
181  CHARACTER(*), INTENT(IN) :: file_header, what
182  INTEGER, INTENT(IN) :: time_step
183  INTEGER :: unit_file=789, j
184  CHARACTER(len=5) :: tit, tit_part
185 
186  IF (what=='new') THEN
187  OPEN (unit=unit_file, file=file_header//'.pvd', form = 'formatted', &
188  access = 'append', status = 'replace')
189  WRITE(unit_file, '(A)') '<?xml version="1.0"?>'
190  WRITE(unit_file, '(A)') '<VTKFile type="Collection" version="0.1"'// &
191  ' byte_order="LittleEndian" compressor="vtkZLibDataCompressor">'
192  WRITE(unit_file, '(A)') '<Collection>'
193  ELSE
194  OPEN (unit=unit_file, file=file_header//'.pvd', form = 'formatted', &
195  access = 'append', status = 'old')
196  backspace(unit_file)
197  backspace(unit_file)
198  END IF
199  WRITE(tit,'(I5)') time_step
200  DO j = 1, SIZE(file_list)
201  WRITE(tit_part,'(I5)') j
202  WRITE(unit_file,'(A)') '<DataSet timestep="'//trim(adjustl(tit))//&
203  '" group="" part="'// &
204  trim(adjustl(tit_part))//'" file="./'//trim(adjustl(file_list(j)))//&
205  '.vtu'//'"/>'
206  END DO
207  WRITE(unit_file, '(A)') '</Collection>'
208  WRITE(unit_file, '(A)') '</VTKFile>'
209 
210  CLOSE(unit_file) !Ecriture pour paraview
211  END SUBROUTINE create_pvd_file
212 
213  SUBROUTINE create_vtu_scal_file(field, mesh, file_name, field_name)
214  USE def_type_mesh
215  IMPLICIT NONE
216  TYPE(mesh_type) :: mesh
217  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: field
218  CHARACTER(*), INTENT(IN) :: file_name, field_name
219  INTEGER :: unit_file=789, m, i, type_cell
220 
221  IF (SIZE(field)==0) RETURN
222 
223  OPEN (unit=unit_file, file=file_name//'.vtu',&
224  form = 'formatted', status = 'unknown')
225 
226  WRITE(unit_file,'(A)') '<VTKFile type="UnstructuredGrid" version="0.1"'// &
227  ' byte_order="LittleEndian">'
228  WRITE(unit_file,'(A)') '<UnstructuredGrid>'
229 
230  WRITE(unit_file,'(A,I9,A,I9,A)') '<Piece NumberOfPoints="', mesh%np, &
231  '" NumberOfCells="', mesh%me, '">'
232  ! PointData Block ------------------------------------------------------------
233  WRITE(unit_file,'(A)') '<PointData Scalars="truc">'
234 !!$ WRITE(unit_file,'(A)') '<DataArray type="Float32" &
235 !!$ Name="'//TRIM(ADJUSTL(field_name))//'" format="ascii">'
236  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="'//trim(adjustl(field_name))//'" format="ascii">'
237  DO i = 1, mesh%np
238  WRITE(unit_file,'(e14.7)') field(i)
239  ENDDO
240  WRITE(unit_file,'(A)') '</DataArray>'
241  WRITE(unit_file,'(A)') '</PointData>'
242  ! End of PointData Block -----------------------------------------------------
243 
244  ! CellData Block -------------------------------------------------------------
245  WRITE(unit_file,'(A)') '<CellData>'
246  WRITE(unit_file,'(A)') '</CellData>'
247  ! End of CellData Block ------------------------------------------------------
248 
249  ! Points Block ---------------------------------------------------------------
250  WRITE(unit_file,'(A)') '<Points>'
251  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="Points" '//&
252  'NumberOfComponents="3" format="ascii">'
253  DO i=1, mesh%np
254  WRITE(unit_file,'(e14.7,A,e14.7)') mesh%rr(1,i), &
255  ' 0.0 ', mesh%rr(2,i)
256  ENDDO
257  WRITE(unit_file,'(A)') '</DataArray>'
258  WRITE(unit_file,'(A)') '</Points>'
259  ! End of Points Block --------------------------------------------------------
260 
261  ! Cells Block ----------------------------------------------------------------
262  IF (mesh%gauss%n_w==3) THEN
263  type_cell = 5
264  ELSE IF (mesh%gauss%n_w==6) THEN
265  type_cell = 22
266  END IF
267  WRITE(unit_file,'(A)') '<Cells>'
268  WRITE(unit_file,'(A)') '<DataArray type="Int64" Name="connectivity" format="ascii">'
269  DO m = 1, mesh%me
270  WRITE(unit_file,'(3(I8,1x))') mesh%jj(1:3,m)-1
271  IF (type_cell==22) THEN
272  WRITE(unit_file,'(3(I8,1x))') mesh%jj(6,m)-1 , mesh%jj(4,m)-1 , mesh%jj(5,m)-1
273  END IF
274  END DO
275  WRITE(unit_file,'(A)') '</DataArray>'
276  WRITE(unit_file,'(A)') '<DataArray type="Int64" Name="offsets" format="ascii">'
277  DO m = 1, mesh%me
278  WRITE(unit_file,'(I8)') m*mesh%gauss%n_w
279  END DO
280  WRITE(unit_file,'(A)') '</DataArray>'
281  WRITE(unit_file,'(A)') '<DataArray type="UInt8" Name="types" format="ascii">'
282  DO m = 1, mesh%me
283  WRITE(unit_file,'(I8)') type_cell
284  END DO
285  WRITE(unit_file,'(A)') '</DataArray>'
286  WRITE(unit_file,'(A)') '</Cells>'
287  ! End of Cells Block ---------------------------------------------------------
288 
289  WRITE(unit_file,'(A)') '</Piece>'
290  WRITE(unit_file,'(A)') '</UnstructuredGrid>'
291  WRITE(unit_file,'(A)') '</VTKFile>'
292 
293  CLOSE(unit_file)
294  END SUBROUTINE create_vtu_scal_file
295 
296  SUBROUTINE create_vtu_vect_file(field, mesh, file_name, opt_st)
297  USE def_type_mesh
298  IMPLICIT NONE
299  TYPE(mesh_type) :: mesh
300  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: field
301  CHARACTER(*), INTENT(IN) :: file_name
302  INTEGER :: unit_file=789, m, i, type_cell
303  CHARACTER(*), OPTIONAL, INTENT(IN) :: opt_st
304  CHARACTER(LEN=200) :: field_name
305 
306  IF (present(opt_st)) THEN
307  field_name=opt_st
308  ELSE
309  field_name='field'
310  END IF
311 
312  IF (SIZE(field)==0) RETURN
313  OPEN (unit=unit_file, file=file_name//'.vtu',&
314  form = 'formatted', status = 'unknown')
315 
316  WRITE(unit_file,'(A)') '<VTKFile type="UnstructuredGrid" version="0.1"'// &
317  ' byte_order="LittleEndian">'
318  WRITE(unit_file,'(A)') '<UnstructuredGrid>'
319 
320  WRITE(unit_file,'(A,I9,A,I9,A)') '<Piece NumberOfPoints="', mesh%np, &
321  '" NumberOfCells="', mesh%me, '">'
322  ! PointData Block ------------------------------------------------------------
323  WRITE(unit_file,'(A)') '<PointData>'
324  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="'//trim(adjustl(field_name))//&
325  '_cos" format="ascii" NumberOfComponents="3">'
326  DO i = 1, mesh%np
327  WRITE(unit_file,'(e14.7)') field(i,1), field(i,3), field(i,5)
328  ENDDO
329  WRITE(unit_file,'(A)') '</DataArray>'
330  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="'//trim(adjustl(field_name))//&
331  '_sin" format="ascii" NumberOfComponents="3">'
332  DO i = 1, mesh%np
333  WRITE(unit_file,'(e14.7)') field(i,2), field(i,4), field(i,6)
334  ENDDO
335  WRITE(unit_file,'(A)') '</DataArray>'
336  WRITE(unit_file,'(A)') '</PointData>'
337  ! End of PointData Block -----------------------------------------------------
338 
339  ! CellData Block -------------------------------------------------------------
340  WRITE(unit_file,'(A)') '<CellData>'
341  WRITE(unit_file,'(A)') '</CellData>'
342  ! End of CellData Block ------------------------------------------------------
343 
344  ! Points Block ---------------------------------------------------------------
345  WRITE(unit_file,'(A)') '<Points>'
346  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="Points" '//&
347  'NumberOfComponents="3" format="ascii">'
348  DO i=1, mesh%np
349  WRITE(unit_file,'(e14.7,A,e14.7)') mesh%rr(1,i), ' 0.0 ' , &
350  mesh%rr(2,i)
351  ENDDO
352  WRITE(unit_file,'(A)') '</DataArray>'
353  WRITE(unit_file,'(A)') '</Points>'
354  ! End of Points Block --------------------------------------------------------
355 
356  ! Cells Block ----------------------------------------------------------------
357  IF (mesh%gauss%n_w==3) THEN
358  type_cell = 5
359  ELSE IF (mesh%gauss%n_w==6) THEN
360  type_cell = 22
361  END IF
362  WRITE(unit_file,'(A)') '<Cells>'
363  WRITE(unit_file,'(A)') '<DataArray type="Int64" Name="connectivity" format="ascii">'
364  DO m = 1, mesh%me
365  WRITE(unit_file,'(3(I8,1x))') mesh%jj(1:3,m)-1
366  IF (type_cell==22) THEN
367  WRITE(unit_file,'(3(I8,1x))') mesh%jj(6,m)-1 , mesh%jj(4,m)-1 , mesh%jj(5,m)-1
368  END IF
369  END DO
370  WRITE(unit_file,'(A)') '</DataArray>'
371  WRITE(unit_file,'(A)') '<DataArray type="Int64" Name="offsets" format="ascii">'
372  DO m = 1, mesh%me
373  WRITE(unit_file,'(I8)') m*mesh%gauss%n_w
374  END DO
375  WRITE(unit_file,'(A)') '</DataArray>'
376  WRITE(unit_file,'(A)') '<DataArray type="UInt8" Name="types" format="ascii">'
377  DO m = 1, mesh%me
378  WRITE(unit_file,'(I8)') type_cell
379  END DO
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_vtu_vect_file
390 
391  SUBROUTINE make_vtu_file_2d(communicator, mesh, header, field, field_name, what, opt_it)
392  USE def_type_mesh
393  USE my_util
394  IMPLICIT NONE
395  TYPE(mesh_type), INTENT(IN) :: mesh
396  CHARACTER(*), INTENT(IN) :: header
397  CHARACTER(*), INTENT(IN) :: field_name, what
398  INTEGER, OPTIONAL, INTENT(IN) :: opt_it
399  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: field
400  INTEGER :: j, it
401  CHARACTER(LEN=200), DIMENSION(:), POINTER :: file_list
402  CHARACTER(LEN=3) :: st_rank, st_it
403 #include "petsc/finclude/petsc.h"
404  petscerrorcode :: ierr
405  petscmpiint :: rank, nb_procs
406  mpi_comm :: communicator
407  CALL mpi_comm_rank(communicator, rank, ierr)
408  CALL mpi_comm_size(communicator, nb_procs, ierr)
409  ALLOCATE(file_list(nb_procs))
410  IF (present(opt_it)) THEN
411  it = opt_it
412  WRITE(st_it,'(I3)') it
413  DO j = 1, nb_procs
414  WRITE(st_rank,'(I3)') j
415  file_list(j) = trim(header)//'_proc_'//trim(adjustl(st_rank))//&
416  '_it_'//trim(adjustl(st_it))
417  END DO
418  ELSE
419  DO j = 1, nb_procs
420  WRITE(st_rank,'(I3)') j
421  file_list(j) = trim(header)//'_proc_'//trim(adjustl(st_rank))
422  END DO
423  END IF
424 
425  CALL check_list(communicator, file_list, mesh%np)
426  IF (rank==0) THEN
427  IF (present(opt_it)) THEN
428  it = opt_it
429  ELSE
430  it = 1
431  END IF
432  CALL create_pvd_file(file_list, trim(header), it, trim(what))
433  END IF
434 !=========TEST FL Feb. 11th, 2013
435  IF (SIZE(field,2) == 6) THEN
436  CALL create_vtu_vect_file(field, mesh, trim(adjustl(file_list(rank+1))), field_name)
437  ELSE IF (SIZE(field,2) == 1) THEN
438  CALL create_vtu_scal_file(field(:,1), mesh, trim(adjustl(file_list(rank+1))), field_name)
439  ELSE IF (SIZE(field,2) .GT. 0) THEN
440  CALL create_vtu_scal_file(field(:,1), mesh, trim(adjustl(file_list(rank+1))), field_name)
441  ELSE
442  CALL error_petsc('Bug in make_vtu_file_2D: field needs at least one component')
443  END IF
444 !=========TEST FL Feb. 11th, 2013
445  END SUBROUTINE make_vtu_file_2d
446 
447  SUBROUTINE make_vtu_file_arpack(communicator, mesh, header, field, field_name, what, num_vp)
448  USE def_type_mesh
449  USE my_util
450  IMPLICIT NONE
451  TYPE(mesh_type), INTENT(IN) :: mesh
452  CHARACTER(*), INTENT(IN) :: header
453  CHARACTER(*), INTENT(IN) :: field_name, what
454  INTEGER, INTENT(IN) :: num_vp
455  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: field
456  CHARACTER(LEN=200), DIMENSION(1) :: file_list
457  CHARACTER(LEN=3) :: st_rank
458 #include "petsc/finclude/petsc.h"
459  petscerrorcode :: ierr
460  petscmpiint :: rank, nb_procs
461  mpi_comm :: communicator
462  CALL mpi_comm_rank(communicator, rank, ierr)
463  CALL mpi_comm_size(communicator, nb_procs, ierr)
464 
465  WRITE(st_rank,'(I3)') num_vp
466  file_list(1) = trim(header)//'_eigen_'//trim(adjustl(st_rank))
467 
468  CALL create_pvd_file(file_list, trim(header), num_vp, trim(what))
469 
470 !=========TEST FL Feb. 11th, 2013
471  IF (SIZE(field,2) == 6) THEN
472  CALL create_vtu_vect_file(field, mesh, trim(adjustl(file_list(1))), field_name)
473  ELSE IF (SIZE(field,2) == 1) THEN
474  CALL create_vtu_scal_file(field(:,1), mesh, trim(adjustl(file_list(1))), field_name)
475  ELSE IF (SIZE(field,2) .GT. 0) THEN
476  CALL create_vtu_scal_file(field(:,1), mesh, trim(adjustl(file_list(1))), field_name)
477  ELSE
478  CALL error_petsc('Bug in make_vtu_file_arpack: field needs at least one component')
479  END IF
480 !=========TEST FL Feb. 11th, 2013
481  END SUBROUTINE make_vtu_file_arpack
482 
483  SUBROUTINE make_vtu_file_axi3d(communicator, mesh, header, &
484  field, field_name, what, opt_it)
485  USE def_type_mesh
486  USE my_util
487  IMPLICIT NONE
488  TYPE(mesh_type), INTENT(IN) :: mesh
489  CHARACTER(*), INTENT(IN) :: header
490  CHARACTER(*), INTENT(IN) :: field_name, what
491  INTEGER, OPTIONAL, INTENT(IN) :: opt_it
492  REAL(KIND=8), DIMENSION(:,:,:),INTENT(IN) :: field
493  INTEGER :: j, it
494  CHARACTER(LEN=200), DIMENSION(:), POINTER :: file_list
495  CHARACTER(LEN=3) :: st_rank
496 #include "petsc/finclude/petsc.h"
497  petscerrorcode :: ierr
498  petscmpiint :: rank, nb_procs
499  mpi_comm :: communicator
500  CALL mpi_comm_rank(communicator, rank, ierr)
501  CALL mpi_comm_size(communicator, nb_procs, ierr)
502  ALLOCATE(file_list(nb_procs))
503  DO j = 1, nb_procs
504  WRITE(st_rank,'(I3)') j
505  file_list(j) = trim(header)//'_proc_'//trim(adjustl(st_rank))
506  END DO
507  CALL check_list(communicator, file_list, mesh%np)
508  IF (rank==0) THEN
509  IF (present(opt_it)) THEN
510  it = opt_it
511  ELSE
512  it = 1
513 
514  END IF
515  CALL create_pvd_file(file_list, trim(header), it, trim(what))
516  END IF
517 
518  CALL create_vtu_file_axi3d(field, mesh, trim(adjustl(file_list(rank+1))), &
519  opt_st=field_name)
520  END SUBROUTINE make_vtu_file_axi3d
521 
522  SUBROUTINE create_vtu_file_axi3d(field, mesh, file_name, opt_st)
523  USE def_type_mesh
524  IMPLICIT NONE
525  TYPE(mesh_type) :: mesh
526  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: field
527  CHARACTER(*), INTENT(IN) :: file_name
528  CHARACTER(*), OPTIONAL, INTENT(IN) :: opt_st
529  CHARACTER(LEN=200) :: field_name
530  REAL(KIND=8) :: theta, pi, dtheta
531  INTEGER :: unit_file=789, m, i, type_cell, &
532  nb_angle, k
533 
534  IF (SIZE(field,2)==0) RETURN
535 
536  IF (present(opt_st)) THEN
537  field_name=opt_st
538  ELSE
539  field_name='field'
540  END IF
541 
542  nb_angle = SIZE(field,1)
543  pi = acos(-1.d0)
544  dtheta = 2*pi/nb_angle
545 
546  OPEN (unit=unit_file, file=file_name//'.vtu',&
547  form = 'formatted', status = 'unknown')
548 
549  WRITE(unit_file,'(A)') '<VTKFile type="UnstructuredGrid" version="0.1"'// &
550  ' byte_order="LittleEndian">'
551  WRITE(unit_file,'(A)') '<UnstructuredGrid>'
552 
553  WRITE(unit_file,'(A,I9,A,I9,A)') '<Piece NumberOfPoints="', nb_angle*mesh%np, &
554  '" NumberOfCells="', nb_angle*mesh%me, '">'
555  ! PointData Block ------------------------------------------------------------
556  WRITE(unit_file,'(A)') '<PointData Scalars="truc">'
557 
558  IF (SIZE(field,2)==1) THEN
559  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="'//trim(adjustl(field_name))//&
560  '" format="ascii">'
561  DO k = 1, nb_angle
562  DO i = 1, mesh%np
563  WRITE(unit_file,'(e14.7)') field(k, 1, i)
564  ENDDO
565  END DO
566  ELSE
567  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="'//trim(adjustl(field_name))//&
568  '" format="ascii" NumberOfComponents="3">'
569  DO k = 1, nb_angle
570  DO i = 1, mesh%np
571  WRITE(unit_file,'(3(e14.7,x))') field(k, 1, i), field(k, 2, i), field(k, 3, i)
572  ENDDO
573  END DO
574  END IF
575 
576  WRITE(unit_file,'(A)') '</DataArray>'
577  WRITE(unit_file,'(A)') '</PointData>'
578  ! End of PointData Block -----------------------------------------------------
579 
580  ! CellData Block -------------------------------------------------------------
581  WRITE(unit_file,'(A)') '<CellData>'
582  WRITE(unit_file,'(A)') '</CellData>'
583  ! End of CellData Block ------------------------------------------------------
584 
585  ! Points Block ---------------------------------------------------------------
586  WRITE(unit_file,'(A)') '<Points>'
587  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="Points" '//&
588  'NumberOfComponents="3" format="ascii">'
589  DO k = 1, nb_angle
590  theta = (k-1)*dtheta
591  DO i=1, mesh%np
592  WRITE(unit_file,'(3(e14.7,2x))') mesh%rr(1,i)*cos(theta), &
593  mesh%rr(1,i)*sin(theta), mesh%rr(2,i)
594  ENDDO
595  END DO
596  WRITE(unit_file,'(A)') '</DataArray>'
597  WRITE(unit_file,'(A)') '</Points>'
598  ! End of Points Block --------------------------------------------------------
599 
600  ! Cells Block ----------------------------------------------------------------
601  type_cell = 13
602  WRITE(unit_file,'(A)') '<Cells>'
603  WRITE(unit_file,'(A)') '<DataArray type="Int64" Name="connectivity" format="ascii">'
604  DO k = 1, nb_angle-1
605  DO m = 1, mesh%me
606  WRITE(unit_file,'(3(I8,1x))') mesh%jj(1:3,m)-1+(k-1)*mesh%np
607  WRITE(unit_file,'(3(I8,1x))') mesh%jj(1:3,m)-1+k*mesh%np
608  END DO
609  END DO
610  k = nb_angle
611  DO m = 1, mesh%me
612  WRITE(unit_file,'(3(I8,1x))') mesh%jj(1:3,m)-1+(k-1)*mesh%np
613  WRITE(unit_file,'(3(I8,1x))') mesh%jj(1:3,m)-1
614  END DO
615  WRITE(unit_file,'(A)') '</DataArray>'
616  WRITE(unit_file,'(A)') '<DataArray type="Int64" Name="offsets" format="ascii">'
617  DO m = 1, nb_angle*mesh%me
618  WRITE(unit_file,'(I8)') 6*m
619  END DO
620  WRITE(unit_file,'(A)') '</DataArray>'
621  WRITE(unit_file,'(A)') '<DataArray type="UInt8" Name="types" format="ascii">'
622  DO m = 1, nb_angle*mesh%me
623  WRITE(unit_file,'(I8)') type_cell
624  END DO
625  WRITE(unit_file,'(A)') '</DataArray>'
626  WRITE(unit_file,'(A)') '</Cells>'
627  ! End of Cells Block ---------------------------------------------------------
628 
629  WRITE(unit_file,'(A)') '</Piece>'
630  WRITE(unit_file,'(A)') '</UnstructuredGrid>'
631  WRITE(unit_file,'(A)') '</VTKFile>'
632 
633  CLOSE(unit_file)
634  END SUBROUTINE create_vtu_file_axi3d
635 
636 
637  SUBROUTINE make_vtu_file_3d(communicator, mesh, header, &
638  field, field_name, what, opt_it)
639  USE def_type_mesh
640  USE my_util
641  IMPLICIT NONE
642  TYPE(mesh_type), INTENT(IN) :: mesh
643  CHARACTER(*), INTENT(IN) :: header
644  CHARACTER(*), INTENT(IN) :: field_name, what
645  INTEGER, OPTIONAL, INTENT(IN) :: opt_it
646  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: field
647  INTEGER :: j, it
648  CHARACTER(LEN=200), DIMENSION(:), POINTER :: file_list
649  CHARACTER(LEN=3) :: st_rank, st_it
650 #include "petsc/finclude/petsc.h"
651  petscerrorcode :: ierr
652  petscmpiint :: rank, nb_procs
653  mpi_comm :: communicator
654  CALL mpi_comm_rank(communicator, rank, ierr)
655  CALL mpi_comm_size(communicator, nb_procs, ierr)
656  ALLOCATE(file_list(nb_procs))
657 
658  IF (present(opt_it)) THEN
659  it = opt_it
660  WRITE(st_it,'(I3)') it
661  DO j = 1, nb_procs
662  WRITE(st_rank,'(I3)') j
663  file_list(j) = trim(header)//'_proc_'//trim(adjustl(st_rank))//&
664  '_it_'//trim(adjustl(st_it))
665  END DO
666  ELSE
667  DO j = 1, nb_procs
668  WRITE(st_rank,'(I3)') j
669  file_list(j) = trim(header)//'_proc_'//trim(adjustl(st_rank))
670  END DO
671  END IF
672 
673  CALL check_list(communicator, file_list, mesh%np)
674  IF (rank==0) THEN
675  IF (present(opt_it)) THEN
676  it = opt_it
677  ELSE
678  it = 1
679 
680  END IF
681  CALL create_pvd_file(file_list, trim(header), it, trim(what))
682  END IF
683 
684  CALL create_vtu_file_3d(field, mesh, trim(adjustl(file_list(rank+1))), &
685  opt_st=field_name)
686  END SUBROUTINE make_vtu_file_3d
687 
688 
689  SUBROUTINE create_vtu_file_3d(field, mesh, file_name, opt_st)
690  USE def_type_mesh
691  USE my_util
692  IMPLICIT NONE
693  TYPE(mesh_type) :: mesh
694  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: field
695  CHARACTER(*), INTENT(IN) :: file_name
696  CHARACTER(*), OPTIONAL, INTENT(IN) :: opt_st
697  CHARACTER(LEN=200) :: field_name
698  INTEGER :: unit_file=789, m, n, type_cell, stride
699 
700  IF (SIZE(field,2)==0) RETURN
701 
702  IF (SIZE(mesh%jj,1)==6) THEN
703  type_cell = 13
704  stride = 6
705  ELSE iF (SIZE(mesh%jj,1)==15) THEN
706  type_cell = 26
707  stride = 15
708  ELSE
709  CALL error_petsc('Bug in create_vtu_file_3D: SIZE(mesh%jj,1) is wrong')
710  END IF
711 
712  IF (present(opt_st)) THEN
713  field_name=opt_st
714  ELSE
715  field_name='field'
716  END IF
717 
718  OPEN (unit=unit_file, file=file_name//'.vtu',&
719  form = 'formatted', status = 'unknown')
720 
721  WRITE(unit_file,'(A)') '<VTKFile type="UnstructuredGrid" version="0.1"'// &
722  ' byte_order="LittleEndian">'
723  WRITE(unit_file,'(A)') '<UnstructuredGrid>'
724 
725  WRITE(unit_file,'(A,I9,A,I9,A)') '<Piece NumberOfPoints="', mesh%np, &
726  '" NumberOfCells="', mesh%me, '">'
727  ! PointData Block ------------------------------------------------------------
728  WRITE(unit_file,'(A)') '<PointData Scalars="truc">'
729 
730  IF (SIZE(field,1)==1) THEN
731  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="'//trim(adjustl(field_name))//&
732  '" format="ascii">'
733  DO n = 1, mesh%np
734  WRITE(unit_file,'(e14.7)') field(1,n)
735  END DO
736  ELSE
737  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="'//trim(adjustl(field_name))//&
738  '" format="ascii" NumberOfComponents="3">'
739  DO n = 1, mesh%np
740  WRITE(unit_file,'(3(e14.7,x))') field(1,n), field(2,n), field(3,n)
741  END DO
742  END IF
743 
744  WRITE(unit_file,'(A)') '</DataArray>'
745  WRITE(unit_file,'(A)') '</PointData>'
746  ! End of PointData Block -----------------------------------------------------
747 
748  ! CellData Block -------------------------------------------------------------
749  WRITE(unit_file,'(A)') '<CellData>'
750  WRITE(unit_file,'(A)') '</CellData>'
751  ! End of CellData Block ------------------------------------------------------
752 
753  ! Points Block ---------------------------------------------------------------
754  WRITE(unit_file,'(A)') '<Points>'
755  WRITE(unit_file,'(A)') '<DataArray type="Float32" Name="Points" '//&
756  'NumberOfComponents="3" format="ascii">'
757  DO n = 1, mesh%np
758  WRITE(unit_file,'(3(e14.7,x))') mesh%rr(1,n), mesh%rr(2,n), mesh%rr(3,n)
759  END DO
760  WRITE(unit_file,'(A)') '</DataArray>'
761  WRITE(unit_file,'(A)') '</Points>'
762  ! End of Points Block --------------------------------------------------------
763 
764  ! Cells Block ----------------------------------------------------------------
765  WRITE(unit_file,'(A)') '<Cells>'
766  WRITE(unit_file,'(A)') '<DataArray type="Int64" Name="connectivity" format="ascii">'
767  DO m = 1, mesh%me
768  WRITE(unit_file,'(15(I8,1x))') mesh%jj(:,m)-1
769  END DO
770  WRITE(unit_file,'(A)') '</DataArray>'
771  WRITE(unit_file,'(A)') '<DataArray type="Int64" Name="offsets" format="ascii">'
772  DO m = 1, mesh%me
773  WRITE(unit_file,'(I8)') stride*m
774  END DO
775  WRITE(unit_file,'(A)') '</DataArray>'
776  WRITE(unit_file,'(A)') '<DataArray type="UInt8" Name="types" format="ascii">'
777  DO m = 1, mesh%me
778  WRITE(unit_file,'(I8)') type_cell
779  END DO
780  WRITE(unit_file,'(A)') '</DataArray>'
781  WRITE(unit_file,'(A)') '</Cells>'
782  ! End of Cells Block ---------------------------------------------------------
783 
784  WRITE(unit_file,'(A)') '</Piece>'
785  WRITE(unit_file,'(A)') '</UnstructuredGrid>'
786  WRITE(unit_file,'(A)') '</VTKFile>'
787 
788  CLOSE(unit_file)
789  END SUBROUTINE create_vtu_file_3d
790 
791  FUNCTION nb_digit() RESULT(dg)
792  USE my_util
793  IMPLICIT NONE
794  INTEGER :: code, nb_procs, dg
795 #include "petsc/finclude/petsc.h"
796  CALL mpi_comm_size(petsc_comm_world,nb_procs,code)
797 
798  IF (nb_procs>9999) THEN
799  CALL error_petsc('nb_procs>9999')
800  END IF
801 
802  IF (nb_procs < 10) THEN
803  dg=1
804  ELSE IF (nb_procs < 99) THEN
805  dg=2
806  ELSE IF (nb_procs < 999) THEN
807  dg=3
808  ELSE IF (nb_procs < 9999) THEN
809  dg=4
810  ELSE
811  dg=0
812  END IF
813  END FUNCTION nb_digit
814 
815 
816  SUBROUTINE make_vtu_file_scalar_2d(comm, mesh, header, field, field_name, what, opt_it)
817  USE def_type_mesh
818  USE my_util
819  IMPLICIT NONE
820  TYPE(mesh_type), INTENT(IN) :: mesh
821  CHARACTER(*), INTENT(IN) :: header
822  CHARACTER(*), INTENT(IN) :: field_name, what
823  INTEGER, OPTIONAL, INTENT(IN) :: opt_it
824  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: field
825  INTEGER :: j, it
826  CHARACTER(LEN=200), DIMENSION(:), POINTER :: file_list
827  CHARACTER(LEN=3) :: st_rank, st_it
828 #include "petsc/finclude/petsc.h"
829  petscerrorcode :: ierr
830  petscmpiint :: rank, nb_procs
831  mpi_comm :: comm
832  CALL mpi_comm_rank(comm, rank, ierr)
833  CALL mpi_comm_size(comm, nb_procs, ierr)
834  ALLOCATE(file_list(nb_procs))
835  IF (present(opt_it)) THEN
836  it = opt_it
837  WRITE(st_it,'(I3)') it
838  DO j = 1, nb_procs
839  WRITE(st_rank,'(I3)') j
840  file_list(j) = trim(header)//'_proc_'//trim(adjustl(st_rank))//&
841  '_it_'//trim(adjustl(st_it))
842  END DO
843  ELSE
844  DO j = 1, nb_procs
845  WRITE(st_rank,'(I3)') j
846  file_list(j) = trim(header)//'_proc_'//trim(adjustl(st_rank))
847  END DO
848  END IF
849 
850  CALL check_list(comm, file_list, mesh%np)
851  IF (rank==0) THEN
852  IF (present(opt_it)) THEN
853  it = opt_it
854  ELSE
855  it = 1
856  END IF
857  CALL create_pvd_file(file_list, trim(header), it, trim(what))
858  END IF
859  CALL create_vtu_scal_file(field, mesh, trim(adjustl(file_list(rank+1))), trim(adjustl(field_name)))
860  END SUBROUTINE make_vtu_file_scalar_2d
861 
862 END MODULE vtk_viz
vtk_viz::create_vtu_file_3d
subroutine, public create_vtu_file_3d(field, mesh, file_name, opt_st)
Definition: plot_vtk.f90:689
vtk_viz::create_vtu_file_axi3d
subroutine, public create_vtu_file_axi3d(field, mesh, file_name, opt_st)
Definition: plot_vtk.f90:522
vtk_viz::create_vtu_vect_file
subroutine, public create_vtu_vect_file(field, mesh, file_name, opt_st)
Definition: plot_vtk.f90:296
chaine_caractere::eval_blank
integer function eval_blank(len_str, string)
Definition: chaine_caractere_par.f90:32
vtk_viz
Definition: plot_vtk.f90:1
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:73
theta
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 theta
Definition: doc_intro.h:193
vtk_viz::make_vtu_file_2d
subroutine, public make_vtu_file_2d(communicator, mesh, header, field, field_name, what, opt_it)
Definition: plot_vtk.f90:391
def_type_mesh
Definition: def_type_mesh.f90:4
my_util
Definition: my_util.f90:1
vtk_viz::make_vtu_file_axi3d
subroutine, public make_vtu_file_axi3d(communicator, mesh, header, field, field_name, what, opt_it)
Definition: plot_vtk.f90:483
vtk_viz::check_list
subroutine, public check_list(communicator, file_list, check)
Definition: plot_vtk.f90:145
vtk_viz::make_vtu_file_scalar_2d
subroutine, public make_vtu_file_scalar_2d(comm, mesh, header, field, field_name, what, opt_it)
Definition: plot_vtk.f90:816
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:447
chaine_caractere
Definition: chaine_caractere_par.f90:4
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:637
vtk_viz::nb_digit
integer function nb_digit()
Definition: plot_vtk.f90:791
vtk_viz::create_pvd_file
subroutine, public create_pvd_file(file_list, file_header, time_step, what)
Definition: plot_vtk.f90:178
my_util::error_petsc
subroutine error_petsc(string)
Definition: my_util.f90:15
vtk_viz::create_vtk_file
subroutine, public create_vtk_file(comm, field, mesh, file_name, opt_it)
Definition: plot_vtk.f90:14
vtk_viz::create_vtu_scal_file
subroutine, public create_vtu_scal_file(field, mesh, file_name, field_name)
Definition: plot_vtk.f90:213
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:193