program termoelastico_1D

!     Versão original 1.0 (21 Mai 02)

!     última alteração =   12 Out 02

!     MODELO NUMÉRICO (resumo)
!       método numérico: diferenças finitas
!       funções de interpolação: CDS, DDS-2 e UDS-2
!       distância entre os nós: constante
!       solver: TDMA
!       precisão: dupla
!       soluções analíticas para L = 1 m

!     ARQUIVOS envolvidos no programa:
!       termoelastico_1D.f90 = programa principal
!       dados.ent            = arquivo de dados do programa
!       saida.txt            = listagem dos resultados
!       t.dat, u., e., s.,   = arquivos de dados para fazer gráfico
!       t.gnu, u., e., s.,   = arquivo de comandos para gerar gráfico
!       notepad.exe          = editor dos arquivos
!       Wgnuplot.exe         = visualizador do gráfico

! -----------------------------------------------

use msfwin

use msflib

implicit none

integer :: n   ! número total de nós

integer :: i   ! número do nó
! i = 1, nó do contorno esquerdo
! i = n, nó do contorno direito
! 2 <= i <= n-1, nós internos

real*8  :: c1, c2 ! auxiliares

real*8  :: dx   ! distância entre nós consecutivos
real*8  :: alfa ! coeficiente de expansão térmica
real*8  :: k    ! condutividade térmica
real*8  :: E    ! módulo de elasticidade
real*8  :: qp   ! geração de calor por unidade de volume
real*8  :: T0   ! temperatura no contorno esquerdo
real*8  :: TL   ! temperatura no contorno direito
real*8  :: Ax   ! área da seção transversal na direção X
real*8  :: L    ! comprimento do domínio de cálculo

real*8  :: F0n  ! solução numérica da força no contorno esquerdo
real*8  :: FLn  ! solução numérica da força no contorno direito

real*8  :: F0a  ! solução análitica da força no contorno esquerdo
real*8  :: FLa  ! solução analítica da força no contorno direito

real*8,dimension(:),allocatable   :: Tn   ! solução numérica da temperatura
real*8,dimension(:),allocatable   :: Ta   ! solução analítica da temperatura

real*8,dimension(:),allocatable   :: un   ! solução numérica do deslocamento
real*8,dimension(:),allocatable   :: ua   ! solução analítica do deslocamento

real*8,dimension(:),allocatable   :: en   ! solução numérica da deformação
real*8,dimension(:),allocatable   :: ea   ! solução analítica da deformação

real*8,dimension(:),allocatable   :: sn   ! solução numérica da tensão
real*8,dimension(:),allocatable   :: sa   ! solução analítica da tensão

real*8,dimension(:),allocatable   :: x    ! coordenada espacial nodal

real*8,dimension(:,:),allocatable :: a    ! coeficiente de f
real*8,dimension(:),  allocatable :: b    ! fonte de f
! primeiro índice = número do ponto
! segundo  índice = conexão com ponto;  1 = j-1;   2 = j;   3 = j+1

integer(4) :: fim ! fecha framewindow

logical :: dos  ! acessa prompt dos
logical :: name ! nome da framewindow

character*100 :: title ! título do gráfico
character*112 :: head  ! título do gráfico + dia

character*12 :: dia       ! data da simulação
character*8  :: hora      ! horário da simulação
integer*4    :: var(8)    ! data e hora
character*20 :: vardate,vartime,varzone ! data e hora
character*70 :: note_caso ! notepad + caso
character*2  :: aux1,aux2
character*4  :: aux3
character*50 :: aux

!-------------------------------------------------

  name=setwindowtext(gethwndqq(qwin$framewindow),  &
       "TERMOELÁSTICO 1D"C)

! -----------------------------------------------

  ! obtém data e hora do sistema

  call date_and_time(vardate,vartime,varzone,var)

  write(aux,*) var(3)
  aux1 = trim(adjustl(aux))
  write(aux,*) var(2)
  aux2 = trim(adjustl(aux))
  write(aux,*) var(1)
  aux3 = trim(adjustl(aux))
  dia = '('//trim(aux1)//'/'//trim(aux2)//'/'//aux3//')'

  write(aux,*) var(5)
  aux1 = trim(adjustl(aux))
  write(aux,*) var(6)
  aux2 = trim(adjustl(aux))
  write(aux,*) var(7)
  aux3 = trim(adjustl(aux))
  hora = trim(aux1)//':'//trim(aux2)//':'//aux3

! -----------------------------------------------

  ! lê dados

  dos = systemqq('notepad dados.ent') ! lista dados

  open(7,file='dados.ent')

    read(7,*) N
    read(7,*) L
    read(7,*) T0
    read(7,*) TL
    read(7,*) k
    read(7,*) qp
    read(7,*) alfa
    read(7,*) E
    read(7,*) Ax
    read(7,*) title

  close(7)

  head = trim(title)//" "//trim(dia)

! -----------------------------------------------

  ! cria o arquivo de saída
  open(10,file='saida.txt')

  write(10,18) head, hora
  18 format(/,5x,'Caso e Dia = ',a112,/, &
              5x,'Hora = ',a8)

  ! alocação de memória
  allocate ( x(n), Tn(n), Ta(n), un(n), ua(n), en(n), ea(n), sn(n), sa(n) )
  allocate ( a(n,3), b(n) )
 
  Tn = 0.0d0
  un = 0.0d0
  en = 0.0d0
  sn = 0.0d0



  ! cálculo da coordenada X

  dx = L / (N-1)
  x(1) = 0.0d0
  do i = 2, N-1
     x(i)  = (i-1) * dx
  end do
  x(N) = 1.0d0



  ! cálculo das soluções analíticas

  c1 = alfa * (TL - T0 + qp/(2*k))
  c2 = - alfa * qp / k  
  do i = 1, N
     Ta(i) = T0 + ( TL - T0 + qp * (1.0d0 - X(i))/(2*k) ) * X(i)
     ua(i) = ( 3.0d0*c1*(X(i)-1.0d0) + c2*((X(i)**2)-1.0d0) ) * X(i) / 6.0d0
     ea(i) = (3.0d0*c1*(2.0d0*x(i)-1.0d0) + c2*(3.0d0*(X(i)**2)-1.0d0))/6.0d0
     sa(i) = E * (ea(i) - alfa*(Ta(i)-T0) )
  end do
  F0a = E * alfa * Ax * ( 3.0d0*(T0-TL) - qp/(2.0d0*k) ) / 6.0d0
  FLa = F0a





  ! cálculo dos coeficientes da Temperatura (T)

    ! contorno esquerdo
    a(1,1) = 0.0d0
    a(1,3) = 0.0d0
    a(1,2) = 1.0d0
    b(1)   = T0

    ! contorno direito
    a(N,1) = 0.0d0
    a(N,3) = 0.0d0
    a(N,2) = 1.0d0
    b(N)   = TL

    ! pontos internos
    do i = 2, N-1
       a(i,1) =  1.0d0
       a(i,3) =  1.0d0
       a(i,2) =  2.0d0
       b(i)   = (qp / k) * (dx**2)
    end do

  ! solução do sistema de equações lineares de T

  call tdma (N,a,b,Tn)

  open(7,file='T.dat')

    write(10,1)
    1 format(/,10x,'TEMPERATURAS (T)', &
             /,t4,'X',t28,'T analítico',t52,'T numérico',/)

    do i = 1, N
       write( 7,2) X(i), Ta(i), Tn(i)
       write(10,2) X(i), Ta(i), Tn(i)
       2 format(3(1pe24.15))
    end do

  close(7)

  open(7,file='T.gnu')
     do k = 1, 6
        read(7,*)
     end do
     write(7,3) head
     3 format("set title '",a112,/, &
              "replot")
  close(7)

  ! plotagem de gráfico com soluções analíticas e numéricas de T

  dos = systemqq('wgnuplot T.gnu')





  ! cálculo dos coeficientes do deslocamento (u)

    ! contorno esquerdo
    a(1,1) = 0.0d0
    a(1,3) = 0.0d0
    a(1,2) = 1.0d0
    b(1)   = 0.0d0

    ! contorno direito
    a(N,1) = 0.0d0
    a(N,3) = 0.0d0
    a(N,2) = 1.0d0
    b(N)   = 0.0d0 

    ! pontos internos
    do i = 2, N-1
       a(i,1) =  1.0d0
       a(i,3) =  1.0d0
       a(i,2) =  2.0d0
       b(i)   = -(alfa * (Tn(i+1)-Tn(i-1)) * dx ) / 2.0d0
    end do

  ! solução do sistema de equações de u

  call tdma (N,a,b,un)

  open(7,file='u.dat')

    write(10,11)
    11 format(/,10x,'DESLOCAMENTOS (u)', &
              /,t4,'X',t28,'u analítico',t52,'u numérico',/)

    do i = 1, N
       write( 7,2) X(i), ua(i), un(i)
       write(10,2) X(i), ua(i), un(i)
    end do

  close(7)

  open(7,file='u.gnu')
     do k = 1, 6
        read(7,*)
     end do
     write(7,3) head
  close(7)

  ! plotagem de gráfico com soluções analíticas e numéricas de u

  dos = systemqq('wgnuplot u.gnu')


 


  ! cálculo das deformações (e)

    ! contorno esquerdo
    en(1) = (4*un(2) - 3*un(1) - un(3))/(2*dx)

    ! contorno direito
    en(N) = (-4*un(N-1) + 3*un(N) + un(N-2))/(2*dx)

    ! pontos internos
    do i = 2, N-1
      en(i) = (un(i+1) - un(i-1))/(2*dx)
    end do

  open(7,file='e.dat')

    write(10,21)
    21 format(/,10x,'DEFORMAÇÕES (e)', &
              /,t4,'X',t28,'e analítico',t52,'e numérico',/)

    do i = 1, N
       write( 7,2) X(i), ea(i), en(i)
       write(10,2) X(i), ea(i), en(i)
    end do

  close(7)

  open(7,file='e.gnu')
     do k = 1, 6
        read(7,*)
     end do
     write(7,3) head
  close(7)

  ! plotagem de gráfico com soluções analíticas e numéricas de e

  dos = systemqq('wgnuplot e.gnu')





  ! cálculo das tensões (s)

    do i = 1, N
      sn(i) = E * (en(i) - alfa*(Tn(i)-Tn(1)))
    end do

  open(7,file='s.dat')

    write(10,31)
    31 format(/,10x,'TENSÕES (s)', &
              /,t4,'X',t28,'s analítico',t52,'s numérico',/)

    do i = 1, N
       write( 7,2) X(i), sa(i), sn(i)
       write(10,2) X(i), sa(i), sn(i)
    end do

  close(7)

  open(7,file='s.gnu')
     do k = 1, 6
        read(7,*)
     end do
     write(7,3) head
  close(7)

  ! plotagem de gráfico com soluções analíticas e numéricas de s

  dos = systemqq('wgnuplot s.gnu')




  F0n = sn(1) * Ax

  FLn = sn(N) * Ax

  write(10,41) F0a, F0n, FLa, FLn
  41 format(/,5x,'FORÇAS NOS CONTORNOS',/,  &
            1pe24.15, ' = esquerdo (analítico)',/, &
            1pe24.15, ' = esquerdo (numérico)',/, &
            1pe24.15, ' = direito (analítico)',/, &
            1pe24.15, ' = direito (numérico)')




  close (10)

  ! abre listagem dos resultados

  note_caso = 'notepad saida.txt'

  dos = systemqq(note_caso) ! lista arquivo de resultados

!  deallocate ( x, fc, fa, a, b )

  fim = setexitqq(qwin$exitnopersist) 

end

!-------------------------------------------------

  subroutine tdma (n,as,bs,xs)

    ! resolve o sistema:  as(i,2)*xs(i) = as(i,1)*xs(i-1) + as(i,3)*xs(i+1) + bs(i)

    implicit none

    integer :: i   ! número do nó
    real*8  :: div ! variável auxiliar

    integer,intent(in) :: n   ! número de nós

    real*8,dimension(:),allocatable :: ps ! coeficiente do tdma
    real*8,dimension(:),allocatable :: qs ! coeficiente do tdma

    real*8,intent(in), dimension(n,3) :: as ! coeficiente do tdma
    real*8,intent(in), dimension(n)   :: bs ! fonte do tdma
    real*8,intent(out),dimension(n)   :: xs ! solução do tdma
    ! primeiro índice = número do volume de controle
    ! segundo  índice = conexão com volume;  1 = W;   2 = P;   3 = E

    allocate(ps(n),qs(n))

    ps(1) = as(1,3) / as(1,2)
    qs(1) = bs(1) / as(1,2)

    do i = 2,n
      div = as(i,2) - as(i,1)*ps(i-1)
      ps(i) = as(i,3) / div
      qs(i) = (bs(i) + as(i,1)*qs(i-1))/div
    end do

    xs(n) = qs(n)

    do i = n-1,1,-1
      xs(i) = ps(i)*xs(i+1) + qs(i)
    end do

    deallocate(ps,qs)

  end subroutine tdma

!-------------------------------------------------

