c==============================================================================
c=====            SUBRUTINA AUXILIAR PARA POSTPROCESO DE DATOS             ====
c==============================================================================
c
c     Esta subrutina genera un fichero de postproceso de resultados en
c     formato *.vtu que puede ser visualizado con ayuda del software
c     Paraview a partir de valores escalares de una determinada magnitud
c     fisica sobre una malla de puntos de discretizacion de un
c     dominio.
c     Los valores del campo escalar a representar se proporcionaran
c     ordenados de izquierda a derecha y de abajo arriba en la malla 2D.
c     Para visualizar el campo escalar en Paraview hay que seleccionarlo
c     en el menu desplegable de la aplicacion.
c     Se pueden generar secuencias de archivos para distintos instantes
c     de tiempo en problemas transitorios que se pueden visualizar
c     conjuntamente cargando simultaneamente todos los ficheros a la vez.
c     De este modo tambien se pueden generar secuencias de imagenes
c     o videos de evolucion de las soluciones obtenidas.
c     Para visualizacion de resultados en 3D en Paraview (adoptando
c     como valor en la tercera dimension el valor escalar proporcionado)
c     se puede utilizar el filtro de datos "Warp By Escalar" de Paraview.
c     Si se estudia solo una porción del dominio debido a posibles simetrias
c     se puede utilizar el filtro de datos "Reflect" de Paraview tantas veces
c     como sea necesario para representar la solución completa del problema.
c
c AUTORES: I. Couceiro, J. Paris, F. Navarrina, I. Colominas
c                                                                 
c        
c FECHA ULTIMA REVISION: 20240412
c
c Copyright: Este programa es propiedad de su autor.
c            No puede utilizarse con fines comerciales sin conocimiento
c            y autorizacion de su autor.
c            Puede utilizarse con fines educativos sin ninguna restriccion.
c
c Definicion de variables de la subrutina:    [I_8 = integer*8, R_8 = real*8]
c                                             [I = input, O = output]
c
c     n  [I_8] [I]   =  Numero de puntos de discretizacion en
c                       horizontal.
c     m  [I_8] [I]   =  Numero de puntos de discretizacion en
c                       vertical.
c     dx [R_8] [I]   =  Separacion entre puntos en direccion horizontal
c     dy [R_8] [I]   =  Separacion entre puntos en direccion vertical
c     u  [R_8] [I]   =  Vector de movimientos para cada uno de los puntos
c                       en los que se discretiza el dominio.
c                       Se almacenan en este vector u por filas.
c     it [I_8] [I]   =  Numero de iteracion para indicar en el nombre del
c                       archivo de Paraview
c
c     NOTA: Se omiten todas las tildes y caracteres especiales en los
c           textos para ajustarse al estandar ASCII y favorecer la
c           compatibilidad del archivo con todos los equipos y sistemas.
c
c===============================================================================
c===============================================================================

      subroutine genera_vtu_iteracion(n,m,dx,dy,u,it)
      implicit real*8(a-h,o-z)
      implicit integer*8(i-n)
      dimension u(n*m)
      character*60 filename
      
      itmax=100000
      itmax=itmax*itmax  ! 10^10
      
      if (n.le.1)then
        write(6,*)'Genera_vtu_tiempos: Numero de puntos en x',
     &            ' insuficiente'
        return
      endif

      if (m.le.1)then
        write(6,*)'Genera_vtu_tiempos: Numero de puntos en y',
     &            ' insuficiente'
        return
      endif

      if (dx.le.0..or.dy.le.0.)then
        write(6,*)'Genera_vtu_tiempos: Las separaciones entre puntos',
     &            ' deben ser positivas'
        return
      endif
      
      if(it.ge.itmax)then
        write(6,*)'Genera_vtu_tiempos: Numero de iteracion',
     &            ' demasiado elevado'
        write(6,*)'No se puede generar el fichero de postproceso',
     &            ' de Paraview'
        return
      endif
      
      write(filename,'(a,i10.10,a)')'Campo_escalar_iter_',it,'.vtu'
      lf=len(filename)
      npoin=n*m
      nelem=(n-1)*(m-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 
        do j=1,m
          do i=1,n
            write(log_vtu,'(3e12.5)')dble(i-1)*dx,dble(j-1)*dy,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,m-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="',
     &                 'Campo_escalar','" format="ascii">'

        do j=1,m
          do i=1,n
            write(log_vtu,'(e12.5)')u((j-1)*n+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
      
