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

      subroutine genera_vtu_triang_escalar(n,alt,u)
      implicit real*8(a-h,o-z)
      implicit integer*8(i-n)
      dimension u(n*(n-1)/2)
      
      if (n.le.1)then
        write(6,*)'Genera_vtu: Numero de puntos',
     &            ' insuficiente'
        call exit(1)
      endif

      if (alt.le.0.)then
        write(6,*)'Genera_vtu_: La altura debe ser positiva'
        call exit(2)
      endif
      

      npoin=n*(n+1)/2
      nelem=(n-1)*n/2
      
      log_vtu=10
      
      open(unit=log_vtu,file='triang_esc_2D.vtu',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">'
! 
        h=alt/dble(n)
        do j=1,n
          do i=1,j
            write(log_vtu,'(3e12.5)')dble(i-1)*h,dble(j-1)*h,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">'
        write(log_vtu,'(4i10)')0,2,1
        do j=3,n
          ia=0
          iup=j*(j-1)/2 ! Faltaria por anadir +1 pero Paraview empieza en 0
          idown=(j-1)*(j-2)/2 ! Faltaria por anadir +1 pero Paraview empieza en 0
          do i=1,j-2
            write(log_vtu,'(4i10)')idown+ia,idown+ia+1,iup+ia+1,iup+ia
            ia=ia+1
          enddo
          write(log_vtu,'(4i10)')idown+ia,iup+ia+1,iup+ia
        enddo
        
        write(log_vtu,'(a)')'        </DataArray>'
!     Offsets
        write(log_vtu,'(a,a)')'<DataArray type="Int64" Name="offsets" ',
     &                        'format="ascii">'
        ipos=0
        do j=2,n
          do i=1,j-2
            ipos=ipos+4
            write(log_vtu,'(i10)')ipos
          enddo
          ipos=ipos+3
          write(log_vtu,'(i10)')ipos
        enddo
        
        write(log_vtu,'(a)')'        </DataArray>'
!    Element type
        write(log_vtu,'(a,a)')'        <DataArray type="UInt8" ',
     &                        'Name="types" format="ascii">'
        do j=2,n
          do i=1,j-2
            write(log_vtu,'(i1)')9   ! Elemento 2D lineal
          enddo
          write(log_vtu,'(i1)')5
        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">'

        write(log_vtu,'(e12.5)')0.d+00
        do j=1,n-1
          do i=1,j
            write(log_vtu,'(e12.5)')u(j*(j-1)/2+i)
          enddo
          write(log_vtu,'(e12.5)')0.d+00
        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