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 regular de puntos de discretizacion de un
c     dominio 1D recto con simetria de revolucion.
c     Se pueden generar secuencias de archivos numerados correlativamente
c     para distintos instantes o soluciones 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 recomienda utilizar el filtro de datos "Warp By Escalar" de Paraview.
c
c AUTORES: I. Couceiro, J. Pari­s, F. Navarrina, I. Colominas
c                                                                 
c        
c FECHA ULTIMA REVISION: 20250317
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 = input, O = output]
c
c  n      [INTEGER*8]    = Numero de puntos totales de discretizacion.
c  ri     [REAL*8]       = Radio interior del dominio axisimetrico
c  re     [REAL*8]       = Radio exterior del dominio axisimetrico
c  u      [REAL*8]       = Vector de valores del parametro escalar a representar
c                          de cada uno de los puntos de discretizacion.
c  nombre [character*50] = Variable tipo character*50 para indicar el nombre
c                          del archivo de Paraview que se va a generar.
c  num    [INTEGER*8]    = Numero de archivo 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_2D_escalar(n,ri,re,u,nombre,num)
      implicit real*8(a-h,o-z)
      implicit integer*8(i-n)
      dimension u(n)
      character*50 nombre
      character*60 filename
      
      
      pi=3.1415926535897932d+00
            
      if (n.le.1)then
        write(6,*)'Genera_vtu_2D_escalar: Numero de puntos',
     &            ' insuficiente'
        return
      endif
      
      if (ri.lt.0.)then
        write(6,*)'Genera_vtu_2D_escalar: El radio interior',
     &            ' debe ser positivo'
        return
      endif
      
      if (re.le.ri)then
        write(6,*)'Genera_vtu_2D_escalar: El radio exerior deber ser',
     &            ' mayor que el interior'
        return
      endif
      
      if (num.lt.1)then
        write(6,*)'Genera_vtu_2D_escalar: El número de archivo',
     &            ' debe ser positivo'
        return
      endif
      
      write(filename,'(2a,i5.5,a)')trim(nombre),'_',num,'.vtu'
      lf=len(filename)
      dx=(re-ri)/dble(n-1)
      nang=61   ! 61 puntos para hacer 60 particiones angulares
      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)
!        dang=2.*pi/dble(nang-1)   ! Descomentar para dibujar la membrana completa
        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)*dx)*coseno,
     &                  (ri+dble(i-1)*dx)*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="',
     &                 'Campo_escalar','" 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
      