Afivo 0.3
All Classes Namespaces Functions Variables Pages
m_vtk.f90
1!> This file is a modification of code found in Lib_VTK_IO
2!>
3!> For the original code, see https://github.com/szaghi/Lib_VTK_IO. I have
4!> extracted the parts that I needed and simplified them a bit. The license for
5!> this file is thus also GPLv3.
6!>
7!> Author Lib_VTK_IO: Stefano Zaghi, modifications: Jannis Teunissen
8!> \todo Document this module
9module m_vtk
10
11 implicit none
12 public
13
14 integer, parameter :: dp = kind(0.0d0)
15 integer, parameter :: sp = kind(0.0)
16 integer, parameter :: buf_len = 200
17 ! Ifort does not yet support new_line("c"), so for now...
18 character, parameter :: endl = char(10)
19 integer, parameter :: indent_spaces = 2
20 integer, parameter :: bytes_i4 = 4
21 integer, parameter :: bytes_r4 = 4
22 integer, parameter :: bytes_r8 = 8
23
24 type vtk_t
25 character(len=100) :: filename
26 logical :: is_open
27 integer :: funit ! xml file
28 integer :: sunit ! scratch file
29 integer :: indent
30 integer :: offset
31 end type vtk_t
32
33contains
34
35 subroutine vtk_ini_xml(vtkf, filename, vtk_type)
36 type(vtk_t), intent(out) :: vtkf
37 character(len=*), intent(in) :: filename, vtk_type
38
39 ! Initialize vtk type
40 open(newunit=vtkf%funit, file=trim(filename), form='UNFORMATTED', &
41 access='STREAM', status='REPLACE')
42 open(newunit=vtkf%sunit, form='UNFORMATTED', &
43 access='STREAM', status='SCRATCH')
44
45 vtkf%filename = filename
46 vtkf%is_open = .true.
47 vtkf%indent = 0
48 vtkf%offset = 0
49
50 write(vtkf%funit) '<?xml version="1.0"?>' // endl
51 write(vtkf%funit) '<VTKFile type="' // trim(vtk_type) // &
52 '" version="0.1" byte_order="LittleEndian">' // endl
53 vtkf%indent = vtkf%indent + indent_spaces
54 end subroutine vtk_ini_xml
55
56 subroutine vtk_unstr_geo_xml(vtkf, coords, n_nodes, n_cells, n_dim, cycl, time)
57 type(vtk_t), intent(inout) :: vtkf
58 real(dp), intent(in) :: coords(:), time
59 real(sp), allocatable :: wr_coords(:)
60 integer, intent(in) :: n_nodes, n_cells, n_dim, cycl
61 character(len=buf_len) :: bufstr
62 integer :: n_bytes, d
63
64 if (n_dim < 1 .or. n_dim > 3) stop "n_dim should be between 1-3"
65
66 ! Always write 3-d point coordinates
67 allocate(wr_coords(3 * n_nodes))
68 wr_coords = 0
69
70 call vtk_dat_xml(vtkf, "FieldData", .true.)
71 write(vtkf%funit) repeat(' ',vtkf%indent) // '<DataArray type="Float64"'// &
72 ' Name="TIME" NumberOfTuples="1" format="ascii">' // endl
73 vtkf%indent = vtkf%indent + indent_spaces
74 write(bufstr, *) repeat(' ',vtkf%indent), time
75 write(vtkf%funit) trim(bufstr) // endl
76 call vtk_dat_xml(vtkf, "DataArray", .false.)
77 write(vtkf%funit) repeat(' ',vtkf%indent) // '<DataArray type="Int32"'// &
78 ' Name="CYCLE" NumberOfTuples="1" format="ascii">' // endl
79 vtkf%indent = vtkf%indent + indent_spaces
80 write(bufstr, *) repeat(' ',vtkf%indent), cycl
81 write(vtkf%funit) trim(bufstr) // endl
82 call vtk_dat_xml(vtkf, "DataArray", .false.)
83 call vtk_dat_xml(vtkf, "FieldData", .false.)
84
85 do d = 1, n_dim
86 wr_coords(d::3) = real(coords(d::n_dim), sp)
87 end do
88
89 write(bufstr, fmt="(A,I0,A,I0,A)") repeat(' ',vtkf%indent) // &
90 '<Piece NumberOfPoints="', n_nodes, '" NumberOfCells="', n_cells, '">'
91 write(vtkf%funit) trim(bufstr) // endl
92 vtkf%indent = vtkf%indent + indent_spaces
93
94 call vtk_dat_xml(vtkf, "Points", .true.)
95
96 write(bufstr, fmt="(A,I0,A,I0,A)") repeat(' ',vtkf%indent) // &
97 '<DataArray type="Float32" NumberOfComponents="3" Name="Points"' // &
98 ' format="appended" offset="', vtkf%offset, '"/>'
99 write(vtkf%funit) trim(bufstr) // endl
100
101 ! Write raw data
102 n_bytes = 3 * n_nodes * bytes_r4
103 vtkf%offset = vtkf%offset + bytes_i4 + n_bytes
104 write(vtkf%sunit) n_bytes, wr_coords
105
106 call vtk_dat_xml(vtkf, "Points", .false.)
107 end subroutine vtk_unstr_geo_xml
108
109 subroutine vtk_unstr_geo_xml_close(vtkf)
110 type(vtk_t), intent(inout) :: vtkf
111 call vtk_dat_xml(vtkf, "Piece", .false.)
112 end subroutine vtk_unstr_geo_xml_close
113
114 subroutine vtk_unstr_con_xml(vtkf, connects, offsets, cell_types, n_cells)
115 type(vtk_t), intent(inout) :: vtkf
116 integer, intent(IN) :: n_cells !< Number of cells.
117 integer, intent(IN) :: connects(:) !< Mesh connectivity.
118 integer, intent(IN) :: offsets(:) !< Cell offset.
119 integer, intent(IN) :: cell_types(:) !< VTK cell type.
120 character(len=buf_len) :: bufstr
121 integer :: n_bytes
122
123 call vtk_dat_xml(vtkf, "Cells", .true.)
124
125 write(bufstr, fmt="(A,I0,A)") repeat(' ',vtkf%indent) // &
126 '<DataArray type="Int32" Name="connectivity" format="appended" offset="', &
127 vtkf%offset, '"/>'
128 write(vtkf%funit) trim(bufstr) // endl
129
130 n_bytes = offsets(n_cells) * bytes_i4
131 vtkf%offset = vtkf%offset + bytes_i4 + n_bytes
132 write(vtkf%sunit) n_bytes, connects
133
134 write(bufstr, fmt="(A,I0,A)") repeat(' ',vtkf%indent) // &
135 '<DataArray type="Int32" Name="offsets" format="appended" offset="', &
136 vtkf%offset, '"/>'
137 write(vtkf%funit) trim(bufstr) // endl
138
139 n_bytes = n_cells * bytes_i4
140 vtkf%offset = vtkf%offset + bytes_i4 + n_bytes
141 write(vtkf%sunit) n_bytes, offsets
142
143 write(bufstr, fmt="(A,I0,A)") repeat(' ',vtkf%indent) // &
144 '<DataArray type="Int32" Name="types" format="appended" offset="', &
145 vtkf%offset, '"/>'
146 write(vtkf%funit) trim(bufstr) // endl
147
148 n_bytes = n_cells * bytes_i4
149 vtkf%offset = vtkf%offset + bytes_i4 + n_bytes
150 write(vtkf%sunit) n_bytes, cell_types
151
152 call vtk_dat_xml(vtkf, "Cells", .false.)
153 end subroutine vtk_unstr_con_xml
154
155 subroutine vtk_dat_xml(vtkf, xml_name, true_is_open)
156 type(vtk_t), intent(inout) :: vtkf
157 character(*), intent(IN) :: xml_name
158 logical, intent(in) :: true_is_open
159
160 if (true_is_open) then
161 write(vtkf%funit) &
162 repeat(' ', vtkf%indent) // '<' // trim(xml_name) // '>' // endl
163 vtkf%indent = vtkf%indent + indent_spaces
164 else
165 vtkf%indent = vtkf%indent - indent_spaces
166 write(vtkf%funit) &
167 repeat(' ', vtkf%indent) // '</' // trim(xml_name) // '>' // endl
168 end if
169 end subroutine vtk_dat_xml
170
171 subroutine vtk_var_r8_xml(vtkf, varname, var, n_data)
172 type(vtk_t), intent(inout) :: vtkf
173 integer, intent(IN) :: n_data !< Number of cells or nodes.
174 character(*), intent(IN) :: varname !< Variable name.
175 real(dp), intent(IN) :: var(:) !< Variable to be saved [1:n_cells_NN].
176 character(len=buf_len) :: bufstr
177 integer :: n_bytes
178
179 write(bufstr, fmt="(A,I0,A)") repeat(' ',vtkf%indent) // &
180 '<DataArray type="Float64" Name="' // trim(varname) // &
181 '" NumberOfComponents="1" format="appended" offset="', &
182 vtkf%offset, '"/>'
183 write(vtkf%funit) trim(bufstr) // endl
184 n_bytes = n_data * bytes_r8
185 vtkf%offset = vtkf%offset + bytes_i4 + n_bytes
186 write(vtkf%sunit) n_bytes, var
187 end subroutine vtk_var_r8_xml
188
189 subroutine vtk_end_xml(vtkf)
190 use, intrinsic :: iso_c_binding
191 type(vtk_t), intent(inout) :: vtkf
192 integer :: n_bytes
193 integer, allocatable :: buffer(:)
194
195 write(vtkf%funit) &
196 repeat(' ',vtkf%indent) // '<AppendedData encoding="raw">' // endl
197 inquire(vtkf%sunit, pos=n_bytes)
198 n_bytes = n_bytes - 1
199 rewind(vtkf%sunit)
200 allocate(buffer(n_bytes/bytes_i4))
201
202 read(vtkf%sunit) buffer(:)
203 write(vtkf%funit) '_', buffer, endl
204 write(vtkf%funit) &
205 repeat(' ',vtkf%indent) // '</AppendedData>' // endl
206 call vtk_dat_xml(vtkf, "VTKFile", .false.)
207 close(vtkf%sunit)
208 close(vtkf%funit)
209 vtkf%is_open = .false.
210 end subroutine vtk_end_xml
211
212end module m_vtk
This file is a modification of code found in Lib_VTK_IO.
Definition: m_vtk.f90:9