Neko  0.9.0
A portable framework for high-order spectral element flow simulations
vtk_file.f90
Go to the documentation of this file.
1 ! Copyright (c) 2019-2022, The Neko Authors
2 ! All rights reserved.
3 !
4 ! Redistribution and use in source and binary forms, with or without
5 ! modification, are permitted provided that the following conditions
6 ! are met:
7 !
8 ! * Redistributions of source code must retain the above copyright
9 ! notice, this list of conditions and the following disclaimer.
10 !
11 ! * Redistributions in binary form must reproduce the above
12 ! copyright notice, this list of conditions and the following
13 ! disclaimer in the documentation and/or other materials provided
14 ! with the distribution.
15 !
16 ! * Neither the name of the authors nor the names of its
17 ! contributors may be used to endorse or promote products derived
18 ! from this software without specific prior written permission.
19 !
20 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21 ! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22 ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23 ! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
24 ! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25 ! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26 ! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
27 ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28 ! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29 ! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30 ! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 ! POSSIBILITY OF SUCH DAMAGE.
32 !
35 module vtk_file
36  use num_types
37  use generic_file
38  use utils
39  use mesh
40  use field
41  use dofmap
42  use mesh_field
43  use tet_mesh
44  use tri_mesh
45  use logger
46  use comm
47  implicit none
48  private
49 
51  type, public, extends(generic_file_t) :: vtk_file_t
52  contains
53  procedure :: read => vtk_file_read
54  procedure :: write => vtk_file_write
55  end type vtk_file_t
56 
57 contains
58 
60  subroutine vtk_file_write(this, data, t)
61  class(vtk_file_t), intent(inout) :: this
62  class(*), target, intent(in) :: data
63  real(kind=rp), intent(in), optional :: t
64  type(mesh_t), pointer :: msh => null()
65  type(field_t), pointer :: fld => null()
66  type(mesh_fld_t), pointer :: mfld => null()
67  type(dofmap_t), pointer :: dm => null()
68  type(tet_mesh_t), pointer :: tet_msh => null()
69  type(tri_mesh_t), pointer :: tri_msh => null()
70  character(len=10) :: id_str
71  integer:: suffix_pos
72 
73  select type(data)
74  type is (mesh_t)
75  msh => data
76  type is(field_t)
77  msh => data%msh
78  fld => data
79  type is(mesh_fld_t)
80  msh => data%msh
81  mfld => data
82  type is (dofmap_t)
83  dm => data
84  type is (tet_mesh_t)
85  tet_msh => data
86  type is (tri_mesh_t)
87  tri_msh => data
88  class default
89  call neko_log%error('Invalid data')
90  end select
91 
92  if (pe_size .gt. 1) then
93  write(id_str,'(i10.10)') pe_rank
94  suffix_pos = filename_suffix_pos(this%fname)
95  open(unit=9, file=trim(this%fname(1:suffix_pos-1))//id_str//'.vtk')
96  else
97  open(unit=9, file=trim(this%fname))
98  end if
99 
100  ! Write legacy header
101  write(9, fmt='(A)') '# vtk DataFile Version 2.0'
102  write(9, fmt='(A)') 'Neko'
103  write(9, fmt='(A)') 'ASCII'
104 
105  if (associated(msh)) then
106  write(9, fmt='(A)') 'DATASET UNSTRUCTURED_GRID'
107 
108  call vtk_file_write_mesh(9, msh)
109 
110  if (associated(mfld)) then
111  call vtk_file_write_cell_data(9, mfld)
112  else if (associated(fld)) then
113  call vtk_file_write_point_data(9, fld)
114  end if
115  else if (associated(dm)) then
116  write(9, fmt='(A)') 'DATASET POLYDATA'
117 
119 
120  call vtk_file_write_dofmap_data(9, dm)
121  else if (associated(tet_msh)) then
122  write(9, fmt='(A)') 'DATASET UNSTRUCTURED_GRID'
123  call vtk_file_write_tet_mesh(9, tet_msh)
124  else if (associated(tri_msh)) then
125  write(9, fmt='(A)') 'DATASET UNSTRUCTURED_GRID'
126  call vtk_file_write_tri_mesh(9, tri_msh)
127  else
128  call neko_error('Invalid data')
129  end if
130 
131  close(9)
132  end subroutine vtk_file_write
133 
134  subroutine vtk_file_read(this, data)
135  class(vtk_file_t) :: this
136  class(*), target, intent(inout) :: data
137 
138  call neko_error('VTK file read not implemented')
139  end subroutine vtk_file_read
140 
142  subroutine vtk_file_write_mesh(unit, msh)
143  integer :: unit
144  type(mesh_t), intent(inout) :: msh
145  integer :: i, j, vtk_type
146  integer, dimension(8), parameter :: vcyc_to_sym = (/1, 2, 4, 3, &
147  5, 6, 8, 7/)
148  ! Dump coordinates
149  write(unit, fmt='(A,I8,A)') 'POINTS', msh%mpts,' double'
150  do i = 1, msh%mpts
151  write(unit, fmt='(F15.8,F15.8,F15.8)') real(msh%points(i)%x,dp)
152  end do
153 
154  ! Dump cells
155  write(unit, fmt='(A,I8,I8)') 'CELLS', msh%nelv, msh%nelv*(msh%npts+1)
156  j = 0
157  do i = 1, msh%nelv
158  write(unit, fmt='(I8,8I8)') msh%npts, &
159  (msh%get_local(msh%elements(i)%e%pts(vcyc_to_sym(j))%p) - 1, &
160  j=1, msh%npts)
161  end do
162 
163  ! Dump cell type for each element
164  write(unit, fmt='(A,I8)') 'CELL_TYPES', msh%nelv
165  vtk_type = 9
166  if (msh%gdim .eq. 3) vtk_type = 12
167  do i = 1, msh%nelv
168  write(unit, fmt='(I2)') vtk_type
169  end do
170 
171  end subroutine vtk_file_write_mesh
172 
174  subroutine vtk_file_write_cell_data(unit, mfld)
175  integer :: unit
176  type(mesh_fld_t), intent(in) :: mfld
177  integer :: i
178 
179  write(unit, fmt='(A,I8)') 'CELL_DATA', mfld%msh%nelv
180  write(unit, fmt='(A,A,A,I8)') 'SCALARS ', trim(mfld%name), ' int', 1
181  write(unit, fmt='(A)') 'LOOKUP_TABLE default'
182 
183  do i = 1, mfld%msh%nelv
184  write(unit, fmt='(I8)') mfld%data(i)
185  end do
186 
187  end subroutine vtk_file_write_cell_data
188 
192  subroutine vtk_file_write_point_data(unit, fld)
193  integer :: unit
194  type(field_t), intent(inout) :: fld
195  real(kind=dp), allocatable :: point_data(:)
196  integer :: i, j, lx, ly, lz, id(8)
197 
198  if ( (fld%Xh%lx - 1 .gt. 1) .or. &
199  (fld%Xh%ly - 1 .gt. 1) .or. &
200  (fld%Xh%lz - 1 .gt. 1)) then
201  call neko_log%warning("Interpolate high-order data onto a low-order mesh")
202  end if
203 
204  write(unit, fmt='(A,I8)') 'POINT_DATA', fld%msh%mpts
205  write(unit, fmt='(A,A,A,I8)') 'SCALARS ', trim(fld%name), ' double', 1
206  write(unit, fmt='(A)') 'LOOKUP_TABLE default'
207 
208  lx = fld%Xh%lx
209  ly = fld%Xh%ly
210  lz = fld%Xh%lz
211  allocate(point_data(fld%msh%mpts))
212 
213  do i = 1, fld%msh%nelv
214  do j = 1, fld%msh%npts
215  id(j) = fld%msh%get_local(fld%msh%elements(i)%e%pts(j)%p)
216  end do
217 
218  point_data(id(1)) = real(fld%x(1,1,1,i),dp)
219  point_data(id(2)) = real(fld%x(lx,1,1,i),dp)
220  point_data(id(3)) = real(fld%x(1,ly,1,i),dp)
221  point_data(id(4)) = real(fld%x(lx,ly,1,i),dp)
222 
223  if (fld%msh%gdim .eq. 3) then
224  point_data(id(5)) = real(fld%x(1,1,lz,i),dp)
225  point_data(id(6)) = real(fld%x(lx,1,lz,i),dp)
226  point_data(id(7)) = real(fld%x(1,ly,lz,i),dp)
227  point_data(id(8)) = real(fld%x(lx,ly,lz,i),dp)
228  end if
229 
230  end do
231 
232  write(unit, *) point_data
233 
234  deallocate(point_data)
235 
236  end subroutine vtk_file_write_point_data
237 
240  integer :: unit
241  type(dofmap_t), intent(inout) :: dm
242  integer :: i,j,k,l
243 
244  write(unit, fmt='(A,I8,A)') 'POINTS', size(dm%x),' double'
245 
246  do i = 1, dm%msh%nelv
247  do l = 1, dm%Xh%lz
248  do k = 1, dm%Xh%ly
249  do j = 1, dm%Xh%lx
250  write(unit, fmt='(F15.8,F15.8,F15.8)') &
251  real(dm%x(j,k,l,i),dp),&
252  real(dm%y(j,k,l,i),dp),&
253  real(dm%z(j,k,l,i),dp)
254  end do
255  end do
256  end do
257  end do
258 
259  write(unit, fmt='(A,I8,I8)') 'VERTICES', size(dm%x), 2*size(dm%x)
260  do i = 1, size(dm%x)
261  write(unit, fmt='(I8,I8)') 1,i-1
262  end do
263 
264 
265  end subroutine vtk_file_write_dofmap_coordinates
266 
268  subroutine vtk_file_write_dofmap_data(unit, dm)
269  integer :: unit
270  type(dofmap_t), intent(inout) :: dm
271  integer :: i, j, k, l
272 
273  write(unit, fmt='(A,I8)') 'POINT_DATA', size(dm%dof)
274  write(unit, fmt='(A,A,A,I8)') 'SCALARS ', 'dof_id', ' integer', 1
275  write(unit, fmt='(A)') 'LOOKUP_TABLE default'
276 
277  do i = 1, dm%msh%nelv
278  do l = 1, dm%Xh%lz
279  do k = 1, dm%Xh%ly
280  do j = 1, dm%Xh%lx
281  write(unit, fmt='(I8)') real(dm%dof(j,k,l,i),dp)
282  end do
283  end do
284  end do
285  end do
286 
287  write(unit, fmt='(A,A,A,I8)') 'SCALARS ', 'shared_dof', ' integer', 1
288  write(unit, fmt='(A)') 'LOOKUP_TABLE default'
289 
290  do i = 1, dm%msh%nelv
291  do l = 1, dm%Xh%lz
292  do k = 1, dm%Xh%ly
293  do j = 1, dm%Xh%lx
294  if (dm%shared_dof(j,k,l,i)) then
295  write(unit, fmt='(I8)') 1
296  else
297  write(unit, fmt='(I8)') 0
298  end if
299  end do
300  end do
301  end do
302  end do
303 
304  end subroutine vtk_file_write_dofmap_data
305 
307  subroutine vtk_file_write_tet_mesh(unit, tet_msh)
308  integer :: unit
309  type(tet_mesh_t), intent(inout) :: tet_msh
310  integer, parameter :: npts = 4
311  integer :: i, j, vtk_type
312 
313  ! Dump coordinates
314  write(unit, fmt='(A,I8,A)') 'POINTS', tet_msh%msh%mpts,' double'
315  do i = 1, tet_msh%msh%mpts
316  write(unit, fmt='(F15.8,F15.8,F15.8)') real(tet_msh%msh%points(i)%x,dp)
317  end do
318 
319  ! Dump cells
320  write(unit, fmt='(A,I8,I8)') 'CELLS', tet_msh%nelv, tet_msh%nelv*(npts+1)
321  j = 0
322  do i = 1, tet_msh%nelv
323  write(unit, fmt='(I8,8I8)') npts, &
324  (tet_msh%msh%get_local(tet_msh%el(i)%pts(j)%p) - 1, &
325  j=1, npts)
326  end do
327 
328  ! Dump cell type for each element
329  write(unit, fmt='(A,I8)') 'CELL_TYPES', tet_msh%nelv
330  vtk_type = 10
331  do i = 1, tet_msh%nelv
332  write(unit, fmt='(I2)') vtk_type
333  end do
334 
335  end subroutine vtk_file_write_tet_mesh
336 
338  subroutine vtk_file_write_tri_mesh(unit, tri_msh)
339  integer :: unit
340  type(tri_mesh_t), intent(inout) :: tri_msh
341  integer, parameter :: npts = 3
342  integer :: i, j, vtk_type
343 
344  ! Dump coordinates
345  write(unit, fmt='(A,I8,A)') 'POINTS', tri_msh%mpts,' double'
346  do i = 1, tri_msh%mpts
347  write(unit, fmt='(F15.8,F15.8,F15.8)') real(tri_msh%points(i)%x,dp)
348  end do
349 
350  ! Dump cells
351  write(unit, fmt='(A,I8,I8)') 'CELLS', tri_msh%nelv, tri_msh%nelv*(npts+1)
352  j = 0
353  do i = 1, tri_msh%nelv
354  write(unit, fmt='(I8,8I8)') npts, &
355  (tri_msh%el(i)%pts(j)%p%id() - 1, j=1, npts)
356  end do
357 
358  ! Dump cell type for each element
359  write(unit, fmt='(A,I8)') 'CELL_TYPES', tri_msh%nelv
360  vtk_type = 5
361  do i = 1, tri_msh%nelv
362  write(unit, fmt='(I2)') vtk_type
363  end do
364 
365  end subroutine vtk_file_write_tri_mesh
366 
367 end module vtk_file
double real
Definition: device_config.h:12
Definition: comm.F90:1
integer pe_rank
MPI rank.
Definition: comm.F90:28
integer pe_size
MPI size of communicator.
Definition: comm.F90:31
Defines a mapping of the degrees of freedom.
Definition: dofmap.f90:35
Defines a field.
Definition: field.f90:34
Module for file I/O operations.
Definition: file.f90:34
Logging routines.
Definition: log.f90:34
type(log_t), public neko_log
Global log stream.
Definition: log.f90:65
Defines a mesh field.
Definition: mesh_field.f90:35
Defines a mesh.
Definition: mesh.f90:34
integer, parameter, public dp
Definition: num_types.f90:9
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
Defines a tetrahedral mesh.
Definition: tet_mesh.f90:35
Defines a triangular surface mesh.
Definition: tri_mesh.f90:35
Utilities.
Definition: utils.f90:35
pure integer function, public filename_suffix_pos(fname)
Find position (in the string) of a filename's suffix.
Definition: utils.f90:56
Legacy VTK file format.
Definition: vtk_file.f90:35
subroutine vtk_file_write(this, data, t)
Write data in legacy VTK.
Definition: vtk_file.f90:61
subroutine vtk_file_write_tet_mesh(unit, tet_msh)
Write a tetrahedral mesh in legacy VTK format.
Definition: vtk_file.f90:308
subroutine vtk_file_read(this, data)
Definition: vtk_file.f90:135
subroutine vtk_file_write_dofmap_data(unit, dm)
Write a dofmap dm data as point data.
Definition: vtk_file.f90:269
subroutine vtk_file_write_tri_mesh(unit, tri_msh)
Write a triangular mesh in legacy VTK format.
Definition: vtk_file.f90:339
subroutine vtk_file_write_dofmap_coordinates(unit, dm)
Write xyz-coordinates of a dofmap dm as points.
Definition: vtk_file.f90:240
subroutine vtk_file_write_cell_data(unit, mfld)
Write a mesh field mfld as cell data.
Definition: vtk_file.f90:175
subroutine vtk_file_write_point_data(unit, fld)
Write a field fld as point data.
Definition: vtk_file.f90:193
subroutine vtk_file_write_mesh(unit, msh)
Write a mesh in legacy VTK format.
Definition: vtk_file.f90:143
A generic file handler.
Interface for legacy VTK files.
Definition: vtk_file.f90:51