c==============================================================================
c=====            AUXILIARY SUBROUTINE FOR DATA POST-PROCESSING            ====
c==============================================================================
c
c    This subroutine generates a post-processing file of results in
c    *.vtu format which can be visualised with the aid of the Paraview
c    software from scalar values of a given physical quantity on a
c    regular grid of discretisation points of a straight 1D domain
c    with revolution symmetry.
c    Sequences of files can be generated for different time instants
c    in transient problems that can be visualised together by
c    simultaneously loading all the files at the same time.
c    In this way, it is also possible to generate sequences of images
c    or videos of the evolution of the solutions obtained.
c    To visualise the results in 3D in Paraview (adopting the scalar
c    value provided as the value in the third dimension) it is
c    recommended to use Paraview's "Warp By Escalar" data filter.

c
c AUTHORS: I. Couceiro, J. Pari­s, F. Navarrina, I. Colominas
c                                                                 
c        
c DATE OF LAST UPDATE: 20241114
c
c Copyright: This programme is the property of its author.
c It may not be used for commercial purposes without the
c knowledge and permission of the author.
c It may be used for educational purposes without restriction.
c
c Variables definition in this subroutine [I = input, O = output]
c
c     n       [I]   =  Number of total discretization points.
c     dr     [R]   =  Distance between radial discretization points
c     ri      [R]   =  Inner radius of the axisymmetric domain
c     u       [R]   =  Vector of values of the scalar parameter to be represented 
c                         for each of the discretisation points.
c     it      [I]   =  Iteration number used to define the name of the
c                      Paraview file
c
c
c===============================================================================
c===============================================================================



      subroutine generate_vtu_timesteps(n,dr,ri,u,it)
      implicit real*8(a-h,o-z)
      implicit integer*8(i-n)
      dimension u(n)
      character*60 filename
      
      pi=3.1415926535897932d+00
      itmax=100000
      itmax=itmax*itmax  ! 10^10
      
      if (n.le.1)then
        write(6,*)'Generate_vtu_timesteps: Insufficient number of',
     &            ' points'
        return
      endif

      if (dr.le.0.)then
        write(6,*)'Generate_vtu_timesteps: The distance between points',
     &            ' must be positive'
        return
      endif
      
      if (ri.le.0.)then
        write(6,*)'Generate_vtu_timesteps: The inner radius',
     &            ' must be positive'
        return
      endif
      
      if(it.ge.itmax)then
        write(6,*)'Generate_vtu_timesteps: Number of iteration',
     &            ' too large'
        write(6,*)'The Paraview postprocessing file',
     &            ' cannot be generated'
        return
      endif
      
      write(filename,'(a,i10.10,a)')'Scalar_field_iter_',it,'.vtu'
      lf=len(filename)
      re=ri+dble(n-1)*dr
      nang=int(pi*(re+ri)/0.02)+1
      npoin=n*nang
      nelem=(n-1)*(nang-1)
      
      log_vtu=10
      
      open(unit=log_vtu,file=filename(1:lf),status='unknown')

        write(log_vtu,'(a,a)')'<VTKFile type="UnstructuredGrid" ',
     &                   'version="0.1" byte_order="LittleEndian">'
        write(log_vtu,'(a)')'  <UnstructuredGrid>'
        write(log_vtu,'(a,i12,a,i12,a)')'    <Piece NumberOfPoints="',
     &                      npoin,'" NumberOfCells="',nelem,'">'
!
!    coordenadas puntuales
!
        write(log_vtu,'(a)')'      <Points>'
        write(log_vtu,'(a,a,i1,a)')'        <DataArray type=',
     &     '"Float64" NumberOfComponents="',3,'" format="ascii">'
c 
        dang=1.25d+00*pi/dble(nang-1)
        do j=1,nang
          coseno=cos(-pi/4.d+00+(j-1)*dang)
          seno=sin(-pi/4.d+00+(j-1)*dang)
          do i=1,n
            write(log_vtu,'(3e12.5)')(ri+dble(i-1)*dr)*coseno,
     &                  (ri+dble(i-1)*dr)*seno,0.d+00
          enddo
        enddo

        write(log_vtu,'(a)')'        </DataArray>'
!
!     Informacion sobre celdas. Conectividad, offset, types
!
        write(log_vtu,'(a)')'      </Points>'
        write(log_vtu,'(a)')'      <Cells>'
        write(log_vtu,'(a)')
     &    '<DataArray type="Int64" Name="connectivity" format="ascii">'
        do j=1,nang-1
          do i=1,n-1
            ipos=(j-1)*n+i-1
            write(log_vtu,'(4i10)')ipos,ipos+1,ipos+n+1,ipos+n
          enddo
        enddo
        
        write(log_vtu,'(a)')'        </DataArray>'
!     Offsets
        write(log_vtu,'(a,a)')'<DataArray type="Int64" Name="offsets" ',
     &                        'format="ascii">'
        do i=1,nelem
          write(log_vtu,'(i10)')4*i
        enddo
        
        write(log_vtu,'(a)')'        </DataArray>'
!    Element type
        write(log_vtu,'(a,a)')'        <DataArray type="UInt8" ',
     &                        'Name="types" format="ascii">'
        do i=1,nelem
          write(log_vtu,'(i1)')9   ! Elemento 2D lineal
        enddo

        write(log_vtu,'(a)')'        </DataArray>'
        write(log_vtu,'(a)')'      </Cells>'
 
!    Scalar field

        write(log_vtu,'(a)')'      <PointData>'
        write(log_vtu,'(5a)')' <DataArray type="Float64" ',
     &                       'NumberOfComponents="1" Name="',
     &                 'Escalar_field','" format="ascii">'

        do j=1,nang
          do i=1,n
            write(log_vtu,'(e12.5)')u(i)
          enddo
        enddo
        write(log_vtu,'(a)')'        </DataArray>'
        write(log_vtu,'(a)')'      </PointData>'

        write(log_vtu,'(a)')'        </Piece>'
        write(log_vtu,'(a)')'  </UnstructuredGrid>'
        write(log_vtu,'(a)')'</VTKFile>'

        close(log_vtu)
      end
      
