From b63c381fcc30fadf8a98006309f4c5726f6d6b07 Mon Sep 17 00:00:00 2001 From: adrianarce-elemwave Date: Wed, 22 Jul 2026 09:20:31 +0000 Subject: [PATCH] FDTD | Feature | Add VTK observation output --- src_output/mapVTKOutput.F90 | 343 +++++++++++++++++++++++++++++++++ src_output/vtkAPI.F90 | 365 ++++++++++++++++++++++++++++++++++++ 2 files changed, 708 insertions(+) create mode 100644 src_output/mapVTKOutput.F90 create mode 100644 src_output/vtkAPI.F90 diff --git a/src_output/mapVTKOutput.F90 b/src_output/mapVTKOutput.F90 new file mode 100644 index 00000000..95d7b0f7 --- /dev/null +++ b/src_output/mapVTKOutput.F90 @@ -0,0 +1,343 @@ +module mapVTKOutput_m + use FDETYPES_m + use outputTypes_m + use outputUtils_m + use directoryUtils_m + use allocationUtils_m + use vtkAPI_m + use volumicProbeUtils_m + use report_m + + implicit none +contains + subroutine init_mapvtk_output(this, lowerBound, upperBound, field, outputTypeExtension, mpidir, problemInfo) + type(mapvtk_output_t), intent(out) :: this + type(cell_coordinate_t), intent(in) :: lowerBound, upperBound + type(problem_info_t), target ,intent(in) :: problemInfo + integer(kind=SINGLE), intent(in) :: mpidir, field + character(len=BUFSIZE), intent(in) :: outputTypeExtension + + integer(kind=SINGLE) :: i + + this%mainCoords = lowerBound + this%auxCoords = upperBound + this%component = field + + this%path = get_output_path() + call store_relevant_coordinates(this, problemInfo) + + contains + + function get_output_path() result(outputPath) + character(len=BUFSIZE) :: probeBoundsExtension, prefixFieldExtension + character(len=BUFSIZE) :: outputPath + probeBoundsExtension = get_coordinates_extension(this%mainCoords, this%auxCoords, mpidir) + prefixFieldExtension = get_prefix_extension(field, mpidir) + outputPath = & + trim(adjustl(outputTypeExtension))//'_'//trim(adjustl(prefixFieldExtension))//'_'//trim(adjustl(probeBoundsExtension)) + return + end function + + subroutine store_media_tag() + end subroutine store_media_tag + end subroutine init_mapvtk_output + + subroutine store_relevant_coordinates(this, problemInfo) + type(mapvtk_output_t), intent(inout) :: this + type(problem_info_t), pointer, intent(in) :: problemInfo + + integer :: i, j, k, field, counter + + counter = 0 + do k = this%mainCoords%Z, this%auxCoords%Z + do j = this%mainCoords%Y, this%auxCoords%Y + do i = this%mainCoords%X, this%auxCoords%X + do field = iEx, iEz + if (isEdge(field, i, j, k, problemInfo)) then + counter = counter + 1 + end if + end do + do field = iHx, iHz + if (isWithinBounds(field, i, j, k, problemInfo)) then + if (isMaterialExceptPML(field, i, j, k, problemInfo)) then + counter = counter + 1 + end if + if (problemInfo%materialTag%getFaceTag(field, i, j, k) < 0 & + .and. (btest(iabs(problemInfo%materialTag%getFaceTag(field, i, j, k)), field - 1)) & + .and. (.not. isPML(field, i, j, k, problemInfo))) then + counter = counter + 1 + end if + end if + end do + end do + end do + end do + + this%nPoints = counter + call alloc_and_init(this%coords, 3, this%nPoints, -99) + call alloc_and_init(this%materialTag, this%nPoints, -1) + call alloc_and_init(this%currentType, this%nPoints, -1) + + counter = 0 + do k = this%mainCoords%Z, this%auxCoords%Z + do j = this%mainCoords%Y, this%auxCoords%Y + do i = this%mainCoords%X, this%auxCoords%X + do field = iEx, iEz + if (isEdge(field, i, j, k, problemInfo)) then + counter = counter + 1 + call writeFaceTagInfo(this, counter, i, j, k, field, problemInfo%materialTag%getFaceTag(field, i, j, k)) + end if + end do + do field = iHx, iHz + if (isWithinBounds(field, i, j, k, problemInfo)) then + if (isMaterialExceptPML(field, i, j, k, problemInfo)) then + counter = counter + 1 + call writeFaceTagInfo(this, counter, i, j, k, field, problemInfo%materialTag%getFaceTag(field, i, j, k)) + end if + if (problemInfo%materialTag%getFaceTag(field, i, j, k) < 0 & + .and. (btest(iabs(problemInfo%materialTag%getFaceTag(field, i, j, k)), field - 1)) & + .and. .not. isPML(field, i, j, k, problemInfo)) then + counter = counter + 1 + call writeFaceTagInfo(this, counter, i, j, k, field, problemInfo%materialTag%getFaceTag(field, i, j, k)) + end if + end if + end do + end do + end do + end do + end subroutine store_relevant_coordinates + + logical function isMaterialExceptPML(field, i, j, k, problemInfo) + integer, intent(in) :: field, i, j, k + type(problem_info_t), intent(in):: problemInfo + isMaterialExceptPML = .not. isMediaVacuum(field, i, j, k, problemInfo) + isMaterialExceptPML = isMaterialExceptPML .and. (.not. isPML(field, i, j, k, problemInfo)) + end function isMaterialExceptPML + + subroutine writeFaceTagInfo(this, counter, i, j, k, field, tag) + type(mapvtk_output_t), intent(inout) :: this + integer, intent(in) :: i, j, k, counter, tag, field + this%coords(1, counter) = i + this%coords(2, counter) = j + this%coords(3, counter) = k + this%currentType(counter) = currentType(field) + this%materialTag(counter) = tag + end subroutine + + subroutine create_geometry_simulation_vtu(this, control, realXGrid, realYGrid, realZGrid) + implicit none + + type(mapvtk_output_t), intent(in) :: this + type(sim_control_t), intent(in) :: control + real(KIND=RKIND), pointer, dimension(:), intent(in) :: realXGrid, realYGrid, realZGrid + + !type(vtk_file) :: vtkOutput + type(vtk_unstructured_grid), target :: ugrid + + integer :: ierr, i, npts, unit + real(RKIND), allocatable :: x(:), y(:), z(:), materialTag(:) + character(len=BUFSIZE) :: info_str + character(len=BUFSIZE) :: metadata_filename, vtuPath + + integer, allocatable :: conn(:), offsets(:), types(:) + integer :: numNodes, numEdges, numQuads + real(kind=RKIND), allocatable :: nodes(:, :) + integer(kind=SINGLE), allocatable :: edges(:, :), quads(:, :) + + call create_folder(this%path, ierr) + vtuPath = join_path(this%path, get_last_component(this%path))//vtuFileExtension + + call createUnstructuredDataForVTU(this%nPoints, this%coords, this%currentType, nodes, edges, quads, numNodes, numEdges, numQuads, control%vtkindex, realXGrid, realYGrid, realZGrid) + call ugrid%add_points(real(nodes, 4)) + + allocate (conn(2*numEdges + 4*numQuads)) + conn(1:2*numEdges) = reshape(edges, [2*numEdges]) + conn(2*numEdges + 1:2*numEdges + 4*numQuads) = reshape(quads, [4*numQuads]) + + allocate (offsets(numEdges + numQuads)) + do i = 1, numEdges + numQuads + if (i <= numEdges) then + if (i == 1) then + offsets(i) = 2 + else + offsets(i) = offsets(i - 1) + 2 + end if + else + if (i == 1) then + offsets(i) = 4 + else + offsets(i) = offsets(i - 1) + 4 + end if + end if + end do + + allocate (types(numEdges + numQuads)) + types(1:numEdges) = 3 + types(numEdges + 1:numEdges + numQuads) = 9 + + call ugrid%add_cell_connectivity(conn, offsets, types) + if (size(offsets) /= numQuads) then + print *, "Problema con offsets" + end if + if (size(types) /= numQuads) then + print *, "Problema con types" + end if + if (offsets(numQuads) /= size(conn)) then + print *, "Tenemos un problema con conn y offset" + end if + + call ugrid%write_file(vtuPath) + + !--------------------------------------------- + ! Metadata: write to .txt file + !--------------------------------------------- + info_str = 'PEC=0, already_YEEadvanced_byconformal=5, NOTOUCHNOUSE=6, '// & + 'WIRE=7, WIRE-COLISION=8, COMPO=3, DISPER=1, DIEL=2, SLOT=4, '// & + 'CONF=5/6, OTHER=-1 (ADD +0.5 for borders)' + + metadata_filename = trim(this%path)//'.txt' + open (newunit=unit, file=metadata_filename, status='replace', action='write', iostat=ierr) + if (ierr /= 0) then + print *, 'Error opening metadata file: ', metadata_filename + else + write (unit, '(A)') trim(info_str) + close (unit) + end if + + end subroutine create_geometry_simulation_vtu + + logical function isEdge(campo, iii, jjj, kkk, problemInfo) + integer(4), intent(in) :: campo, iii, jjj, kkk + type(problem_info_t), pointer, intent(in) :: problemInfo + + type(MediaData_t), pointer, dimension(:) :: mData + type(limit_t), pointer, dimension(:) :: problemDimension + + integer(4) :: imed, imed1, imed2, imed3, imed4, contaborde + + mData => problemInfo%materialList + problemDimension => problemInfo%problemDimension + isEdge = .false. + contaborde = 0 + + call get_media_from_coord_and_h_neighbours(campo, iii, jjj, kkk, problemInfo%geometryToMaterialData, imed, imed1, imed2, imed3, imed4) + + if (imed /= 1) then + + if (mData(imed)%is%SGBC) then + + if (mData(imed1)%is%SGBC) then + if (trim(adjustl(mData(imed)%Multiport(1)%MultiportFileZ11)) /= & + trim(adjustl(mData(imed1)%Multiport(1)%MultiportFileZ11))) contaborde = contaborde + 1 + elseif (imed1 /= 1) then + contaborde = contaborde + 1 + end if + + if (mData(imed2)%is%SGBC) then + if (trim(adjustl(mData(imed)%Multiport(1)%MultiportFileZ11)) /= & + trim(adjustl(mData(imed2)%Multiport(1)%MultiportFileZ11))) contaborde = contaborde + 1 + elseif (imed2 /= 1) then + contaborde = contaborde + 1 + end if + + if (mData(imed3)%is%SGBC) then + if (trim(adjustl(mData(imed)%Multiport(1)%MultiportFileZ11)) /= & + trim(adjustl(mData(imed3)%Multiport(1)%MultiportFileZ11))) contaborde = contaborde + 1 + elseif (imed3 /= 1) then + contaborde = contaborde + 1 + end if + + if (mData(imed4)%is%SGBC) then + if (trim(adjustl(mData(imed)%Multiport(1)%MultiportFileZ11)) /= & + trim(adjustl(mData(imed4)%Multiport(1)%MultiportFileZ11))) contaborde = contaborde + 1 + elseif (imed4 /= 1) then + contaborde = contaborde + 1 + end if + + elseif (mData(imed)%is%Multiport) then + + if (mData(imed1)%is%Multiport) then + if (trim(adjustl(mData(imed)%Multiport(1)%MultiportFileZ11)) /= & + trim(adjustl(mData(imed1)%Multiport(1)%MultiportFileZ11))) contaborde = contaborde + 1 + elseif (imed1 /= 1) then + contaborde = contaborde + 1 + end if + + if (mData(imed2)%is%Multiport) then + if (trim(adjustl(mData(imed)%Multiport(1)%MultiportFileZ11)) /= & + trim(adjustl(mData(imed2)%Multiport(1)%MultiportFileZ11))) contaborde = contaborde + 1 + elseif (imed2 /= 1) then + contaborde = contaborde + 1 + end if + + if (mData(imed3)%is%Multiport) then + if (trim(adjustl(mData(imed)%Multiport(1)%MultiportFileZ11)) /= & + trim(adjustl(mData(imed3)%Multiport(1)%MultiportFileZ11))) contaborde = contaborde + 1 + elseif (imed3 /= 1) then + contaborde = contaborde + 1 + end if + + if (mData(imed4)%is%Multiport) then + if (trim(adjustl(mData(imed)%Multiport(1)%MultiportFileZ11)) /= & + trim(adjustl(mData(imed4)%Multiport(1)%MultiportFileZ11))) contaborde = contaborde + 1 + elseif (imed4 /= 1) then + contaborde = contaborde + 1 + end if + + elseif (mData(imed)%is%AnisMultiport) then + + if (mData(imed1)%is%AnisMultiport) then + if (trim(adjustl(mData(imed)%AnisMultiport(1)%MultiportFileZ11)) /= & + trim(adjustl(mData(imed1)%AnisMultiport(1)%MultiportFileZ11))) contaborde = contaborde + 1 + elseif (imed1 /= 1) then + contaborde = contaborde + 1 + end if + + if (mData(imed2)%is%AnisMultiport) then + if (trim(adjustl(mData(imed)%AnisMultiport(1)%MultiportFileZ11)) /= & + trim(adjustl(mData(imed2)%AnisMultiport(1)%MultiportFileZ11))) contaborde = contaborde + 1 + elseif (imed2 /= 1) then + contaborde = contaborde + 1 + end if + + if (mData(imed3)%is%AnisMultiport) then + if (trim(adjustl(mData(imed)%AnisMultiport(1)%MultiportFileZ11)) /= & + trim(adjustl(mData(imed3)%AnisMultiport(1)%MultiportFileZ11))) contaborde = contaborde + 1 + elseif (imed3 /= 1) then + contaborde = contaborde + 1 + end if + + if (mData(imed4)%is%AnisMultiport) then + if (trim(adjustl(mData(imed)%AnisMultiport(1)%MultiportFileZ11)) /= & + trim(adjustl(mData(imed4)%AnisMultiport(1)%MultiportFileZ11))) contaborde = contaborde + 1 + elseif (imed4 /= 1) then + contaborde = contaborde + 1 + end if + + else + if ((imed /= imed1) .and. (imed1 /= 1)) contaborde = contaborde + 1 + if ((imed /= imed2) .and. (imed2 /= 1)) contaborde = contaborde + 1 + if ((imed /= imed3) .and. (imed3 /= 1)) contaborde = contaborde + 1 + if ((imed /= imed4) .and. (imed4 /= 1)) contaborde = contaborde + 1 + end if + + if ((imed1 == 1 .and. imed2 == 1 .and. imed3 == 1 .and. imed4 /= 1) .or. & + (imed2 == 1 .and. imed3 == 1 .and. imed4 == 1 .and. imed1 /= 1) .or. & + (imed3 == 1 .and. imed4 == 1 .and. imed1 == 1 .and. imed2 /= 1) .or. & + (imed4 == 1 .and. imed1 == 1 .and. imed2 == 1 .and. imed3 /= 1) .or. & + (imed1 == 1 .and. imed2 == 1 .and. imed3 == 1 .and. imed4 == 1) .or. & + (contaborde > 0)) isEdge = .true. + + if ((iii > problemDimension(campo)%XE) .or. (jjj > problemDimension(campo)%YE) .or. & + (kkk > problemDimension(campo)%ZE)) isEdge = .false. + + if ((iii < problemDimension(campo)%XI) .or. (jjj < problemDimension(campo)%YI) .or. & + (kkk < problemDimension(campo)%ZI)) isEdge = .false. + + else + isEdge = .false. + end if + + end function + +end module mapVTKOutput_m diff --git a/src_output/vtkAPI.F90 b/src_output/vtkAPI.F90 new file mode 100644 index 00000000..bca99f0b --- /dev/null +++ b/src_output/vtkAPI.F90 @@ -0,0 +1,365 @@ +module vtkAPI_m + implicit none + private + public :: vtk_data_array, vtk_grid, vtk_structured_grid, vtk_unstructured_grid + + !========================== + ! Data array type + !========================== + type :: vtk_data_array + character(len=:), allocatable :: name + character(len=:), allocatable :: type + integer :: num_components + real, allocatable :: data(:) + end type vtk_data_array + + !========================== + ! Base abstract class + !========================== + type, abstract :: vtk_grid + real, allocatable :: points(:, :) + type(vtk_data_array), allocatable :: scalars(:) + type(vtk_data_array), allocatable :: vectors(:) + type(vtk_data_array), allocatable :: cell_scalars(:) + type(vtk_data_array), allocatable :: cell_vectors(:) + contains + procedure(add_points_generic), deferred :: add_points + procedure(add_scalar_generic), deferred :: add_scalar + procedure(add_vector_generic), deferred :: add_vector + procedure(add_cell_scalar_generic), deferred :: add_cell_scalar + procedure(add_cell_vector_generic), deferred :: add_cell_vector + procedure(write_file_generic), deferred :: write_file + end type vtk_grid + + !========================== + ! Structured Grid (VTS) + !========================== + type, extends(vtk_grid) :: vtk_structured_grid + integer :: nx = 0, ny = 0, nz = 0 + contains + procedure :: add_points => add_points_structured + procedure :: add_scalar => add_scalar_structured + procedure :: add_vector => add_vector_structured + procedure :: add_cell_scalar => add_cell_scalar_structured + procedure :: add_cell_vector => add_cell_vector_structured + procedure :: write_file => write_vts_file + end type vtk_structured_grid + + !========================== + ! Unstructured Grid (VTU) + !========================== + type, extends(vtk_grid) :: vtk_unstructured_grid + integer :: num_points = 0 + integer :: num_cells = 0 + integer, allocatable :: connectivity(:) + integer, allocatable :: offsets(:) + integer, allocatable :: types(:) + contains + procedure :: add_points => add_points_unstructured + procedure :: add_scalar => add_scalar_unstructured + procedure :: add_vector => add_vector_unstructured + procedure :: add_cell_scalar => add_cell_scalar_unstructured + procedure :: add_cell_vector => add_cell_vector_unstructured + procedure :: add_cell_connectivity + procedure :: write_file => write_vtu_file + end type vtk_unstructured_grid + + !========================== + ! Generic deferred interfaces + !========================== + abstract interface + subroutine add_points_generic(this, pts) + import :: vtk_grid + class(vtk_grid), intent(inout) :: this + real, intent(in) :: pts(:, :) + end subroutine add_points_generic + + subroutine add_scalar_generic(this, name, data) + import :: vtk_grid + class(vtk_grid), intent(inout) :: this + character(len=*), intent(in) :: name + real, intent(in) :: data(:) + end subroutine add_scalar_generic + + subroutine add_vector_generic(this, name, data) + import :: vtk_grid + class(vtk_grid), intent(inout) :: this + character(len=*), intent(in) :: name + real, intent(in) :: data(:) + end subroutine add_vector_generic + + subroutine add_cell_scalar_generic(this, name, data) + import :: vtk_grid + class(vtk_grid), intent(inout) :: this + character(len=*), intent(in) :: name + real, intent(in) :: data(:) + end subroutine add_cell_scalar_generic + + subroutine add_cell_vector_generic(this, name, data) + import :: vtk_grid + class(vtk_grid), intent(inout) :: this + character(len=*), intent(in) :: name + real, intent(in) :: data(:) + end subroutine add_cell_vector_generic + + subroutine write_file_generic(this, filename) + import :: vtk_grid + class(vtk_grid), intent(in) :: this + character(len=*), intent(in) :: filename + end subroutine write_file_generic + end interface + +contains + !========================== + !==== Structured Grid Methods ==== + !========================== + + subroutine add_points_structured(this, pts) + class(vtk_structured_grid), intent(inout) :: this + real, intent(in) :: pts(:, :) + integer :: npts + npts = this%nx*this%ny*this%nz + if (size(pts, 1) /= 3) error stop 'add_points_structured: first dim must be 3' + if (size(pts, 2) /= npts) error stop 'add_points_structured: wrong number of points' + if (allocated(this%points)) deallocate (this%points) + allocate (this%points(3, npts)) + this%points = pts + end subroutine add_points_structured + + subroutine add_scalar_structured(this, name, data) + class(vtk_structured_grid), intent(inout) :: this + character(len=*), intent(in) :: name + real, intent(in) :: data(:) + call add_array_generic(this%scalars, name, data, 1) + end subroutine add_scalar_structured + + subroutine add_vector_structured(this, name, data) + class(vtk_structured_grid), intent(inout) :: this + character(len=*), intent(in) :: name + real, intent(in) :: data(:) + call add_array_generic(this%vectors, name, data, 3) + end subroutine add_vector_structured + + subroutine add_cell_scalar_structured(this, name, data) + class(vtk_structured_grid), intent(inout) :: this + character(len=*), intent(in) :: name + real, intent(in) :: data(:) + call add_array_generic(this%cell_scalars, name, data, 1) + end subroutine add_cell_scalar_structured + + subroutine add_cell_vector_structured(this, name, data) + class(vtk_structured_grid), intent(inout) :: this + character(len=*), intent(in) :: name + real, intent(in) :: data(:) + call add_array_generic(this%cell_vectors, name, data, 3) + end subroutine add_cell_vector_structured + + !========================== + ! Write VTS + !========================== + subroutine write_vts_file(this, filename) + class(vtk_structured_grid), intent(in) :: this + character(len=*), intent(in) :: filename + integer :: iunit + open (newunit=iunit, file=filename, status='replace', action='write', form='formatted') + write (iunit, *) '' + write (iunit, *) ' ' + write (iunit, *) ' ' + call write_pointdata(iunit, this%scalars, this%vectors) + call write_celldata(iunit, this%cell_scalars, this%cell_vectors) + call write_points(iunit, this%points) + write (iunit, *) ' ' + write (iunit, *) ' ' + write (iunit, *) '' + close (iunit) + end subroutine write_vts_file + + !========================== + !==== Unstructured Grid Methods ==== + !========================== + subroutine add_points_unstructured(this, pts) + real, intent(in) :: pts(:, :) + class(vtk_unstructured_grid), intent(inout) :: this + if (size(pts, 1) /= 3) error stop 'add_points_unstructured: first dim must be 3' + this%num_points = size(pts, 2) + if (allocated(this%points)) deallocate (this%points) + allocate (this%points(3, this%num_points)) + this%points = pts + end subroutine add_points_unstructured + + subroutine add_scalar_unstructured(this, name, data) + character(len=*), intent(in) :: name + real, intent(in) :: data(:) + class(vtk_unstructured_grid), intent(inout) :: this + call add_array_generic(this%scalars, name, data, 1) + end subroutine add_scalar_unstructured + + subroutine add_vector_unstructured(this, name, data) + character(len=*), intent(in) :: name + real, intent(in) :: data(:) + class(vtk_unstructured_grid), intent(inout) :: this + call add_array_generic(this%vectors, name, data, 3) + end subroutine add_vector_unstructured + + subroutine add_cell_scalar_unstructured(this, name, data) + character(len=*), intent(in) :: name + real, intent(in) :: data(:) + class(vtk_unstructured_grid), intent(inout) :: this + call add_array_generic(this%cell_scalars, name, data, 1) + end subroutine add_cell_scalar_unstructured + + subroutine add_cell_vector_unstructured(this, name, data) + character(len=*), intent(in) :: name + real, intent(in) :: data(:) + class(vtk_unstructured_grid), intent(inout) :: this + call add_array_generic(this%cell_vectors, name, data, 3) + end subroutine add_cell_vector_unstructured + + subroutine add_cell_connectivity(this, conn, offsets, types) + class(vtk_unstructured_grid), intent(inout) :: this + integer, intent(in) :: conn(:), offsets(:), types(:) + this%num_cells = size(offsets) + if (allocated(this%connectivity)) deallocate (this%connectivity) + if (allocated(this%offsets)) deallocate (this%offsets) + if (allocated(this%types)) deallocate (this%types) + allocate (this%connectivity(size(conn))) + allocate (this%offsets(size(offsets))) + allocate (this%types(size(types))) + this%connectivity = conn + this%offsets = offsets + this%types = types + end subroutine add_cell_connectivity + + !========================== + ! Write VTU + !========================== + subroutine write_vtu_file(this, filename) + character(len=*), intent(in) :: filename + class(vtk_unstructured_grid), intent(in) :: this + integer :: iunit + open (newunit=iunit, file=filename, status='replace', action='write', form='formatted') + write (iunit, *) '' + write (iunit, *) ' ' + write (iunit, *) ' ' + call write_pointdata(iunit, this%scalars, this%vectors) + call write_celldata(iunit, this%cell_scalars, this%cell_vectors) + call write_cells(iunit, this%connectivity, this%offsets, this%types) + call write_points(iunit, this%points) + write (iunit, *) ' ' + write (iunit, *) ' ' + write (iunit, *) '' + close (iunit) + end subroutine write_vtu_file + +!========================== +!==== Shared helper routines ==== +!========================== + subroutine add_array_generic(array_list, name, data, ncomp) + type(vtk_data_array), allocatable, intent(inout) :: array_list(:) + character(len=*), intent(in) :: name + real, intent(in) :: data(:) + integer, intent(in) :: ncomp + type(vtk_data_array), allocatable :: tmp(:) + integer :: n + if (.not. allocated(array_list)) then + allocate (array_list(1)) + n = 1 + else + n = size(array_list) + 1 + allocate (tmp(n)) + tmp(1:n - 1) = array_list + call move_alloc(tmp, array_list) + end if + array_list(n)%name = name + array_list(n)%type = 'Float32' + array_list(n)%num_components = ncomp + array_list(n)%data = data + end subroutine add_array_generic + + subroutine write_pointdata(iunit, scalars, vectors) + integer, intent(in) :: iunit + type(vtk_data_array), allocatable, intent(in) :: scalars(:) + type(vtk_data_array), allocatable, intent(in) :: vectors(:) + integer :: i + character(len=1024) :: tag + tag = ' 0) tag = trim(tag)//' Scalars="'//trim(scalars(1)%name)//'"' + end if + if (allocated(vectors)) then + if (size(vectors) > 0) tag = trim(tag)//' Vectors="'//trim(vectors(1)%name)//'"' + end if + tag = trim(tag)//'>' + write (iunit, '(A)') trim(tag) + if (allocated(scalars)) then + do i = 1, size(scalars) + write (iunit, '(A)') ' ' + write (iunit, '(1000(F12.6,1X))') scalars(i)%data + write (iunit, '(A)') ' ' + end do + end if + if (allocated(vectors)) then + do i = 1, size(vectors) + write (iunit, '(A)') ' ' + write (iunit, '(1000(F12.6,1X))') vectors(i)%data + write (iunit, '(A)') ' ' + end do + end if + write (iunit, *) ' ' + end subroutine write_pointdata + + subroutine write_celldata(iunit, scalars, vectors) + integer, intent(in) :: iunit + type(vtk_data_array), allocatable, intent(in) :: scalars(:) + type(vtk_data_array), allocatable, intent(in) :: vectors(:) + integer :: i + if (.not. allocated(scalars) .and. .not. allocated(vectors)) then + write (iunit, *) ' ' + write (iunit, *) ' ' + return + end if + write (iunit, *) ' ' + if (allocated(scalars)) then + do i = 1, size(scalars) + write (iunit, '(A)') ' ' + write (iunit, '(1000(F12.6,1X))') scalars(i)%data + write (iunit, '(A)') ' ' + end do + end if + if (allocated(vectors)) then + do i = 1, size(vectors) + write (iunit, '(A)') ' ' + write (iunit, '(1000(F12.6,1X))') vectors(i)%data + write (iunit, '(A)') ' ' + end do + end if + write (iunit, *) ' ' + end subroutine write_celldata + + subroutine write_points(iunit, pts) + integer, intent(in) :: iunit + real, intent(in) :: pts(:, :) + write (iunit, *) ' ' + write (iunit, *) ' ' + write (iunit, '(1000(F12.6,1X))') pts + write (iunit, *) ' ' + write (iunit, *) ' ' + end subroutine write_points + + subroutine write_cells(iunit, conn, offsets, types) + integer, intent(in) :: iunit + integer, intent(in) :: conn(:), offsets(:), types(:) + write (iunit, *) ' ' + write (iunit, *) ' ' + write (iunit, '(1000(I8,1X))') conn + write (iunit, *) ' ' + write (iunit, *) ' ' + write (iunit, '(1000(I8,1X))') offsets + write (iunit, *) ' ' + write (iunit, *) ' ' + write (iunit, '(1000(I3,1X))') types + write (iunit, *) ' ' + write (iunit, *) ' ' + end subroutine write_cells + +end module vtkAPI_m