!======================================================================!
!                                                                      !
!    Software Name : FrontCOMP_cure   Ver. 3.1                         !
!                                                                      !
!      Module Name : m_heat_output                                     !
!      Category    : Heat Analysis                                     !
!                                                                      !
!      Developed based on "FrontSTR" of RSS21 project                  !
!                                                                      !
!                     Written by Yasuji Fukahori,    2006/06/01        !
!                                Tomotaka Ogasawara, 2013/03/26        !
!                                                                      !
!     Contact address :  IIS,The University of Tokyo, CISS             !
!                                                                      !
!    "Composite Material Strength & Reliability Evaluation Simulator"  !
!                                                                      !
!======================================================================!

module m_heat_output

use m_heat_LIB_DFLUX

   contains

!C***
!C***  OUTPUT FILES
!C***
!C----------------------------------------------------------------------*
   subroutine heat_output(hecMESH,hecMAT,fstrHEAT,fstrPARAM)
!C----------------------------------------------------------------------*

      use m_fstr

      implicit none

      type (hecmwST_local_mesh) :: hecMESH
      type (hecmwST_matrix    ) :: hecMAT
      type (fstr_heat         ) :: fstrHEAT
      type (fstr_param        ) :: fstrPARAM

! Local variable
      integer(kind=kint) :: i,inod

!C
!C-- TEMPERATURE
!C
      if(fstrHEAT%iout_list(1) .eq. 1) then
        write(ILOG,*)
        write(ILOG,*) '     Node   Temperature   '
        write(ILOG,*) '--------------------------'
        do i= 1, hecMESH%nn_internal
          inod=fstrPARAM%global_local_id(1,i)
          write (ILOG,'(i8,f12.4)') inod, fstrHEAT%TEMP(i)
        enddo
      end if

!C
!C-- CURE
!C

!!      call heat_nodal_cure3d( hecMESH, fstrHEAT)
!!      call hecmw_update_m_R(hecMESH,fstrHEAT%CURE3D_nod,hecMESH%n_node,1)

      if(fstrHEAT%iout_list(2) .eq. 1) then
        write(ILOG,*)
        write(ILOG,*) '     Node   Cure          '
        write(ILOG,*) '--------------------------'
        do i= 1, hecMESH%nn_internal
          inod=fstrPARAM%global_local_id(1,i)
          write (ILOG,'(i8,f12.4)') inod, fstrHEAT%CURE3D_nod(i)
        enddo
      end if

   end subroutine heat_output

!C----------------------------------------------------------------------*
   subroutine heat_output_sta(hecMESH,hecMAT,fstrHEAT,fstrPARAM    &
                             ,heat_sta_IW1,heat_sta_IW2)
!C----------------------------------------------------------------------*

      use m_fstr

      implicit none

      type (hecmwST_local_mesh) :: hecMESH
      type (hecmwST_matrix    ) :: hecMAT
      type (fstr_heat         ) :: fstrHEAT
      type (fstr_param        ) :: fstrPARAM

! Local variable
      integer(kind=kint) :: i,j,inod,ielem,g1,g2,g3
      integer(kind=kint) :: heat_sta_IW1, heat_sta_IW2 
      real(kind=kreal)   :: bb(27)
!C
!C-- TEMPERATURE
!C
      if(fstrHEAT%iout_list(1) .eq. 1) then
        do i= 1, hecMESH%nn_internal
          inod=fstrPARAM%global_local_id(1,i)
!         write (heat_sta_IW1,'(i10,a,1pe13.4e3)') inod,' , ',fstrHEAT%TEMP(i)
          write (heat_sta_IW1) inod,fstrHEAT%TEMP(i)
        enddo
      end if
!!
      if(fstrHEAT%iout_list(1) .eq. -999) then
        do i= 1, hecMAT%NP    
          write (heat_sta_IW1,'(i10,a,1pe13.4e3)') i,' , ',fstrHEAT%TEMP(i)
        enddo
      end if
!C
!C-- CURE
!C
      if(fstrHEAT%iout_list(2) .eq. 1) then
        do i= 1, hecMESH%nn_internal
          inod=fstrPARAM%global_local_id(1,i)
!         write (heat_sta_IW2,'(i10,a,1pe13.4e3)') inod,' , ',fstrHEAT%CURE3D_nod(i)
          write (heat_sta_IW2) inod,fstrHEAT%CURE3D_nod(i)
        enddo
      end if
!!
      if(fstrHEAT%iout_list(2) .eq. -999) then
        do i= 1, hecMAT%NP    
          write (heat_sta_IW2,'(i10,a,1pe13.4e3)') i,' , ',fstrHEAT%CURE3D_nod(i)
        enddo
      end if
   end subroutine heat_output_sta

!C----------------------------------------------------------------------*
      subroutine heat_nodal_cure3d( hecMESH, fstrHEAT)
!C----------------------------------------------------------------------*
!C
!C NODAL CURE3D
!C
      use m_fstr

      implicit none

      type (hecmwST_local_mesh) :: hecMESH
      type (fstr_heat         ) :: fstrHEAT

!** Local variables
      real(kind=kreal) xx(20),yy(20),zz(20)
      real(kind=kreal) disp(60)
      real(kind=kreal) tt(20),tt0(20)
      integer(kind=kint) nodLOCAL(20)
      real(kind=kreal) array(13,20)

      integer(kind=kint) ig0,ig,iS0,iE0,ik,in,ic,itype,iS,iE,ic_type,nn,icel,j 
      integer(kind=kint) isect,k,i 
      integer(kind=kint) ierror 
      real(kind=kreal) val,ee,pp,rho,alpha,thick,s11,s22,s33,s12,s23,s13,ps,smises 
!** Array for nodal recovery
      real    (kind=KREAL), dimension(:,:), allocatable :: arrayTotal
      integer (kind=KINT), dimension(:), allocatable :: number
      real    (kind=KREAL), dimension(:), allocatable :: temp
!*Allocate array
      allocate ( arrayTotal( hecMESH%n_node, 1  )         ,STAT=ierror) 
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <heat_nodal_cure3d, arrayTotal>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      allocate ( number(hecMESH%n_node)                    ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <heat_nodal_cure3d, number>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
!*ZERO CLEAR 
      arrayTotal=0
      number=0
!C
!C +-------------------------------+
!C | according to ELEMENT TYPE     |
!C +-------------------------------+
      ic=0
      do itype= 1, hecMESH%n_elem_type
        iS= hecMESH%elem_type_index(itype-1) + 1
        iE= hecMESH%elem_type_index(itype  )
        ic_type= hecMESH%elem_type_item(itype)
        if (.not. hecmw_is_etype_solid(ic_type)) cycle
!C** Set number of nodes	
        nn = hecmw_get_max_node(ic_type)
!C element loop
        do icel= iS, iE
!C** node ID
          iS= hecMESH%elem_node_index(icel-1)
          do j=1,nn
            nodLOCAL(j)= hecMESH%elem_node_item (iS+j)
!C** nodal coordinate
            xx(j)=hecMESH%node(3*nodLOCAL(j)-2)
            yy(j)=hecMESH%node(3*nodLOCAL(j)-1)
            zz(j)=hecMESH%node(3*nodLOCAL(j)  )
          enddo  
!C** section  ID
          isect= hecMESH%section_ID(icel)
!C** Create local stiffness
          ic=ic+1
          if (ic_type.eq.341) then
          else if ( ic_type.eq.351) then
          else if ( ic_type.eq.361) then
          else if ( ic_type.eq.342) then
!!!!!            call NRCV_342CURE(ILOG,xx,yy,zz,      &
!!!!!                         fstrHEAT%CURE3D(:,:,:,ic),array)
          else if ( ic_type.eq.352) then
          else if ( ic_type.eq.362) then
          endif
          do j=1,nn
            do k=1,1
              arrayTotal(nodLOCAL(j),k)= arrayTotal(nodLOCAL(j),k)+array(k,j)
            enddo
            number( nodLOCAL(j) )=number( nodLOCAL(j) )+1
          enddo
        enddo
      enddo
!** Average over nodes
      do i=1,hecMESH%n_node
        do j=1,1
          arrayTotal(i,j)=arrayTotal(i,j)/FLOAT(number(i))
        enddo
      enddo
!** Show
!C     do i=1,hecMESH%n_node
!C       j=hecMESH%global_node_ID(i)
!C        write(*,'(i8)') j
!C        write(*,'(''E '',6e11.3)') (arrayTotal(i,k),k=1,6)
!C        write(*,'(''S '',7e11.3)') (arrayTotal(i,k),k=7,13)
!C     enddo
!** Set Array 
      do i=1,hecMESH%n_node
        do j=1,1
          fstrHEAT%CURE3D_nod(1*(i-1)+j)=arrayTotal(i,j)
        enddo
      enddo
!*Deallocate array
      deallocate( arrayTotal                 ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <heat_nodal_cure3d, arrayTotal>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate( number                     ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <heat_nodal_cure3d, number>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      return
      end subroutine heat_nodal_cure3d

!C----------------------------------------------------------------------*
      subroutine heat_post(hecMESH,hecMAT,fstrHEAT)
!C----------------------------------------------------------------------*
      use m_fstr

      implicit none

      type (hecmwST_local_mesh) :: hecMESH
      type (hecmwST_matrix)     :: hecMAT
      type (fstr_heat         ) :: fstrHEAT

!** Local variables
!C Max Min
      real(kind=kreal) Tmax(1),Tmin(1),Amax(1),Amin(1)
      integer(kind=kint) I_Tmax(1),I_Tmin(1),I_Amax(1),I_Amin(1)

      integer(kind=kint) j,i,k,ii,MDOF,ID_area,nd,id,ndof 

!C*-------- solver control -----------*
      logical :: ds = .false. !using Direct Solver or not

! in case of direct solver
      if (hecMAT%Iarray(99) .eq. 2) then
        ds = .true.
      end if

!C Clear
      Tmax=0.0d0
      Tmin=0.0d0
      Amax=0.0d0
      Amin=0.0d0
      I_Tmax=0
      I_Tmin=0
      I_Amax=0
      I_Amin=0

!C*** Temp. & Cure 
      write(ILOG,*) '#### TEMPERATURE & CURE @POST'

      do i= 1, hecMESH%nn_internal
        j=hecMESH%global_node_ID(i)
          do k=1,1 
           if(fstrHEAT%idx_mat_node(1*i-1+k).eq.2) then
            Amax(k)=fstrHEAT%CURE3D_nod(1*i-1+k)
            Amin(k)=fstrHEAT%CURE3D_nod(1*i-1+k)
            I_Amax(k)=j
            I_Amin(k)=j
           end if
          enddo
      enddo
!
      do i= 1, hecMESH%nn_internal
        j=hecMESH%global_node_ID(i)
!C Max & Min
        if( i.eq.1) then
          do k=1,1 
            Tmax(k)=fstrHEAT%TEMP(1*i-1+k)
            Tmin(k)=fstrHEAT%TEMP(1*i-1+k)
            I_Tmax(k)=j
            I_Tmin(k)=j
          enddo
        else
          do k=1,1 
            if( fstrHEAT%TEMP(1*i-1+k) .gt. Tmax(k) ) then
              Tmax(k)=fstrHEAT%TEMP(1*i-1+k)
              I_Tmax(k)=j
            endif
            if( fstrHEAT%TEMP(1*i-1+k) .lt. Tmin(k) ) then
              Tmin(k)=fstrHEAT%TEMP(1*i-1+k)
              I_Tmin(k)=j
            endif
          enddo
          do k=1,1 
           if(fstrHEAT%idx_mat_node(1*i-1+k).eq.2) then
            if( fstrHEAT%CURE3D_nod(1*i-1+k) .gt. Amax(k) ) then
              Amax(k)=fstrHEAT%CURE3D_nod(1*i-1+k)
              I_Amax(k)=j
            endif
            if( fstrHEAT%CURE3D_nod(1*i-1+k) .lt. Amin(k) ) then
              Amin(k)=fstrHEAT%CURE3D_nod(1*i-1+k)
              I_Amin(k)=j
            endif
           end if
          enddo
        endif
      enddo
!!----------------------------------------------------------------------!
      write(ILOG,*) '##### Local Summary @Node :Max/IdMax/Min/IdMin####'
      write(ILOG,*) '//T  ',Tmax(1),I_Tmax(1),Tmin(1),I_Tmin(1)
      write(ILOG,*) '//a  ',Amax(1),I_Amax(1),Amin(1),I_Amin(1)

!C*** Show Summary
!!---- oga2009 --    if (ds .ne. .true.) then !In case of Direct Solver prevent MPI
     if (ds .neqv. .true.) then !In case of Direct Solver prevent MPI
       call hecmw_allREDUCE_R(hecMESH,Tmax,1,hecmw_max)
       call hecmw_allREDUCE_R(hecMESH,Tmin,1,hecmw_min)
       call hecmw_allREDUCE_R(hecMESH,Amax,1,hecmw_max)
       call hecmw_allREDUCE_R(hecMESH,Amin,1,hecmw_min)
     end if 

     if( hecMESH%my_rank .eq. 0 ) then
        write(ILOG,*) '##### Global Summary @Node :Max/Min####'
        write(ILOG,*) '//T  ',Tmax(1),Tmin(1)
        write(ILOG,*) '//a  ',Amax(1),Amin(1)
     endif
!C
      return
      end subroutine heat_post

!CC-------------------------------------------------------------------------------
!CC    cure of cfrp for visualizer
!CC-------------------------------------------------------------------------------
   subroutine cure_output_cfrp ( hecMESH,hecMAT,fstrHEAT )

      use m_fstr

      implicit none
      integer(kind=kint) itype,iS,iE,ic_type,icel,isect,IMAT,in0,nn,i,nodLOCAL,ip,inod
      integer(kind=kint) ic,jp,jnod,isU,ieU,ik,isL,ieL  
      real(kind=kreal)   XX,YY,ZZ,TT,T0,CUR_nod_I,ASECT,S0,SS,THICK


      type (fstr_heat         ) :: fstrHEAT
      type (hecmwST_matrix    ) :: hecMAT
      type (hecmwST_local_mesh) :: hecMESH

      dimension nodLOCAL(20),XX(20),YY(20),ZZ(20),TT(20),T0(20),S0(20),CUR_nod_I(20),SS(400)
!!

      do i= 1, hecMESH%n_node
        fstrHEAT%CURE3D_nod_cfrp(i) = fstrHEAT%CURE3D_nod(i)
      end do
!

!C +-------------------------------+
!C | ELEMENT-by-ELEMENT ASSEMBLING | 
!C | according to ELEMENT TYPE     |
!C +-------------------------------+

      do itype= 1, hecMESH%n_elem_type
        iS= hecMESH%elem_type_index(itype-1) + 1
        iE= hecMESH%elem_type_index(itype  )
        ic_type= hecMESH%elem_type_item(itype)
   
        do icel = iS, iE
          isect = hecMESH%section_ID(icel)
          IMAT = hecMESH%section%sect_mat_ID_item(isect)
          in0 = hecMESH%elem_node_index(icel-1)

          nn = hecmw_get_max_node(ic_type)
          do i = 1, nn
            nodLOCAL(i) = hecMESH%elem_node_item(in0+i)
          enddo

          if( IMAT .eq. 2) then
            do i = 1, nn
              fstrHEAT%idx_mat_node   (nodLOCAL(i)) = 2
            end do
          else if( IMAT .ne. 2) then
            do i = 1, nn
              fstrHEAT%CURE3D_nod_cfrp(nodLOCAL(i)) = 1.0D+00
            end do
          end if

        enddo
      enddo

   end subroutine cure_output_cfrp 

end module m_heat_output
