subroutine output_overall_vtk_fields(self, file_counts, domain, porosity, temperature, si, pressure, water_flux)
implicit none
class(type_output_overall), intent(inout) :: self
integer(int32), intent(in) :: file_counts
type(type_domain), intent(in) :: domain
real(real64), intent(in), optional :: porosity(:)
real(real64), intent(in), optional :: temperature(:)
real(real64), intent(in), optional :: si(:)
real(real64), intent(in), optional :: pressure(:)
type(type_dp_3d), intent(in), optional :: water_flux
integer(int32) :: status
integer(int32) :: unit_num
integer(int32) :: iN, iE, idx, i
character(256) :: output_name
real(real64), allocatable :: original(:), original_vector(:, :)
! Initialize VTK file
write (output_name, self%format_output) trim(self%dir_output_field), "Out_", file_counts, self%file_extension
open (newunit=unit_num, file=output_name, status='replace', action='write', iostat=status)
if (status /= 0) call error_message(931)
write (unit_num, '(a)') "# vtk DataFile Version 2.0"
write (unit_num, '(a)') "Analysis ASCII VTK file"
write (unit_num, '(a)') "ASCII"
write (unit_num, '(a)') "DATASET UNSTRUCTURED_GRID"
write (unit_num, '(a,i0,a)') "POINTS ", self%vtk%num_points, " double"
do iN = 1, self%vtk%num_points
write (unit_num, '(3(es22.15,x))') self%vtk%coordinate%x(iN), self%vtk%coordinate%y(iN), self%vtk%coordinate%z(iN)
end do
write (unit_num, '(a)') ""
write (unit_num, '(a,i0,x,i0,a)') "CELLS ", self%vtk%num_cells, sum(self%vtk%offsets(:)) + self%vtk%num_cells
idx = 1
do iE = 1, self%vtk%num_cells
write (unit_num, '(i0,'//to_string(self%vtk%offsets(iE))//'(x,i0))') self%vtk%offsets(iE), &
self%vtk%connectivities(idx:idx + self%vtk%offsets(iE) - 1)
idx = idx + self%vtk%offsets(iE)
end do
write (unit_num, '(a)') ""
write (unit_num, '(a,i0)') "CELL_TYPES ", self%vtk%num_cells
write (unit_num, '(i0)') self%vtk%cell_types(:)
write (unit_num, '(a)') ""
do i = 1, size(self%variable_names)
if (i == 1) write (unit_num, '(a, i0)') "POINT_DATA ", self%vtk%num_points
select case (self%variable_names(i))
case ("temperature")
if (present(temperature)) then
call allocate_array(original, self%vtk%num_points)
call domain%reordering%to_original_value(temperature, original)
call write_field(unit_num, "Temperature", original)
call deallocate_array(original)
end if
case ("ice_saturation")
if (present(si)) then
call allocate_array(original, self%vtk%num_points)
call domain%reordering%to_original_value(si, original)
call write_field(unit_num, "Si", original)
call deallocate_array(original)
end if
case ("thermal_conductivity")
print *, "Warning: 'thermal_conductivity' is not implemented in VTK output."
case ("volumetric_heat_capacity")
print *, "Warning: 'volumetric_heat_capacity' is not implemented in VTK output."
case ("pressure")
if (present(pressure)) then
call allocate_array(original, self%vtk%num_points)
call domain%reordering%to_original_value(pressure, original)
call write_field(unit_num, "Pressure", original)
call deallocate_array(original)
end if
case ("water_flux")
if (present(water_flux)) then
call allocate_array(original_vector, 3_int32, self%vtk%num_points)
call domain%reordering%to_original_value(water_flux%x, original_vector(:, 1))
call domain%reordering%to_original_value(water_flux%y, original_vector(:, 2))
call domain%reordering%to_original_value(water_flux%z, original_vector(:, 3))
call write_field(unit_num, "waterFlux", original_vector(:, 1), original_vector(:, 2), original_vector(:, 3))
call deallocate_array(original_vector)
end if
case ("hydraulic_conductivity")
print *, "Warning: 'hydraulic_conductivity' is not implemented in VTK output."
end select
end do
end subroutine output_overall_vtk_fields