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

module m_heat_mat_ass_bc_DFLUX
   contains
!C
!C***
!C*** MAT_ASS_DFLUX
!C***
!C
   subroutine heat_mat_ass_bc_DFLUX( hecMESH, hecMAT, fstrHEAT, CTIME, DTIME, BETA, incr, iterALL )

      use m_fstr
      use m_heat_get_amplitude
      use m_heat_LIB_DFLUX
!!    use m_heat_LIB_CAPACITY

      implicit none
      integer(kind=kint) k,icel,jj,ic_type,isect,isuf,iamp,nn,iS,j
      integer(kind=kint) IMAT
      integer(kind=kint) incr,iterALL
      real(kind=kreal)   CTIME,QQ,val,asect,thick
      real(kind=kreal)   DTIME,BETA,ALFA
      type (fstr_heat         ) :: fstrHEAT
      type (hecmwST_matrix    ) :: hecMAT
      type (hecmwST_local_mesh) :: hecMESH
      real(kind=kreal) xx(20),yy(20),zz(20),tt(20),tt0(20),CUR_nod(20),CUR_nod_I(20)
      real(kind=kreal) vect(20),vect0(20),d_tt(20) 
      integer(kind=kint) nodLocal(20)
!C
      ALFA = 1.0 - BETA
!
        do k = 1, fstrHEAT%Q_SUF_tot 

          icel    = fstrHEAT%Q_SUF_elem(k)
           jj = hecMESH%global_elem_ID(icel)
          ic_type = hecMESH%elem_type(icel)
          isect   = hecMESH%section_ID(icel)
          IMAT    = hecMESH%section%sect_mat_ID_item(isect)
          isuf    = fstrHEAT%Q_SUF_surf(k)
          iamp    = fstrHEAT%Q_SUF_ampl(k)

          call heat_get_amplitude( fstrHEAT,iamp,CTIME,QQ )
          val     = fstrHEAT%Q_SUF_val (k) * QQ
!C**
          nn = hecmw_get_max_node(ic_type)
!C** 
          iS = hecMESH%elem_node_index(icel-1)
          do j = 1, nn
            nodLOCAL(j) = hecMESH%elem_node_item(iS+j)
            xx(j) = hecMESH%node( 3*nodLOCAL(j)-2 )
            yy(j) = hecMESH%node( 3*nodLOCAL(j)-1 )
            zz(j) = hecMESH%node( 3*nodLOCAL(j)   )
!!!!!    20111014   -----
            tt0 (j) = fstrHEAT%  TEMP0( nodLOCAL(j) )
            tt  (j) = fstrHEAT%  TEMP ( nodLOCAL(j) )
            if ( iterALL .eq. 1 ) then
              d_tt(j) = fstrHEAT%D_TEMP ( nodLOCAL(j) )
            else if ( iterALL .gt. 1 ) then
              d_tt(j) = fstrHEAT%  TEMP ( nodLOCAL(j) ) - fstrHEAT%TEMP0( nodLOCAL(j) )
            end if

            CUR_nod_I(j) = fstrHEAT%CURE3D_nod_I( nodLOCAL(j)   )
            CUR_nod  (j) = fstrHEAT%CURE3D_nod  ( nodLOCAL(j)   )
!!!!!    20111014   -----
          enddo  
!C**		
          if    ( ic_type.eq.111 ) then
            is = hecMesh%section%sect_R_index(isect)
            asect = hecMESH%section%sect_R_item(is)
            call heat_DFLUX_111(nn,xx,yy,zz,asect,isuf,val,vect)

          elseif( ic_type.eq.231 ) then
            is = hecMesh%section%sect_R_index(isect)
            thick = hecMESH%section%sect_R_item(is)
            call heat_DFLUX_231(nn,xx,yy,zz,thick,isuf,val,vect)

          elseif( ic_type.eq.232 ) then
            is = hecMesh%section%sect_R_index(isect)
            thick = hecMESH%section%sect_R_item(is)
            call heat_DFLUX_232(nn,xx,yy,zz,thick,isuf,val,vect)

          elseif( ic_type.eq.241 ) then
            is = hecMesh%section%sect_R_index(isect)
            thick = hecMESH%section%sect_R_item(is)
            call heat_DFLUX_241(nn,xx,yy,zz,thick,isuf,val,vect)

          elseif( ic_type.eq.242 ) then
            is = hecMesh%section%sect_R_index(isect)
            thick = hecMESH%section%sect_R_item(is)
            call heat_DFLUX_242(nn,xx,yy,zz,thick,isuf,val,vect)

          elseif( ic_type.eq.341 ) then
            call heat_DFLUX_341(nn,xx,yy,zz,isuf,val,vect)

          elseif( ic_type.eq.342 ) then
!!          call heat_DFLUX_342(nn,xx,yy,zz,isuf,val,vect)  !--- oga2009
            call heat_DFLUX_342(nn,xx,yy,zz,isuf,val,vect0                            &
                               ,tt0,IMAT,DTIME,icel,jj                                &
                               ,fstrHEAT%RHOtab(IMAT)    ,fstrHEAT%RHOtemp(IMAT,:)    &
                               ,fstrHEAT%RHOfuncA(IMAT,:),fstrHEAT%RHOfuncB(IMAT,:)   &
                               ,fstrHEAT%idx_cure          &
                               ,fstrHEAT%ccc_rl            &
                               ,CUR_nod_I                  &
                               ,IDBG,incr,iterALL,d_tt)
            call heat_DFLUX_342(nn,xx,yy,zz,isuf,val,vect                             &
                               ,tt,IMAT,DTIME,icel,jj                                 &
                               ,fstrHEAT%RHOtab(IMAT)    ,fstrHEAT%RHOtemp(IMAT,:)    &
                               ,fstrHEAT%RHOfuncA(IMAT,:),fstrHEAT%RHOfuncB(IMAT,:)   &
                               ,fstrHEAT%idx_cure          &
                               ,fstrHEAT%ccc_rl            &
                               ,CUR_nod                    &
                               ,IDBG,incr,iterALL,d_tt)

          elseif( ic_type.eq.351 ) then
            call heat_DFLUX_351(nn,xx,yy,zz,isuf,val,vect)

          elseif( ic_type.eq.352 ) then
            call heat_DFLUX_352(nn,xx,yy,zz,isuf,val,vect)

          elseif( ic_type.eq.361 ) then
            call heat_DFLUX_361(nn,xx,yy,zz,isuf,val,vect)

          elseif( ic_type.eq.362 ) then
            call heat_DFLUX_362(nn,xx,yy,zz,isuf,val,vect)

          elseif( ic_type.eq.731 ) then
            is    = hecMesh%section%sect_R_index(isect)
            thick = hecMESH%section%sect_R_item(is)
            call heat_DFLUX_731(nn,xx,yy,zz,thick,isuf,val,vect)

          elseif( ic_type.eq.741 ) then
            is    = hecMesh%section%sect_R_index(isect)
            thick = hecMESH%section%sect_R_item(is)
            call heat_DFLUX_741(nn,xx,yy,zz,thick,isuf,val,vect)

          endif

!C 
          do j = 1, nn
!           hecMAT%B( nodLOCAL(j) ) = hecMAT%B( nodLOCAL(j) ) - vect(j) 
            hecMAT%B( nodLOCAL(j) ) = hecMAT%B( nodLOCAL(j) ) - vect(j)*BETA - vect0(j)*ALFA 
          enddo
!C

        enddo

      return
   end subroutine heat_mat_ass_bc_DFLUX 

!C----------------------------------------------------------------
!C*** MAT_ASS_DFLUX FOR CURE
!C----------------------------------------------------------------
   subroutine cure_mat_ass_bc_DFLUX( hecMESH, hecMAT, fstrHEAT, CTIME, DTIME, BETA, incr,iterCUR )

      use m_fstr
      use m_heat_get_amplitude
      use m_heat_LIB_DFLUX
!!    use m_heat_LIB_CAPACITY

      implicit none
      integer(kind=kint) k,icel,jj,ic_type,isect,isuf,iamp,nn,iS,j
      integer(kind=kint) IMAT
      integer(kind=kint) incr,iterCUR
      real(kind=kreal)   CTIME,QQ,val,asect,thick
      real(kind=kreal)   DTIME,BETA,ALFA
      type (fstr_heat         ) :: fstrHEAT
      type (hecmwST_matrix    ) :: hecMAT
      type (hecmwST_local_mesh) :: hecMESH
      real(kind=kreal) xx(20),yy(20),zz(20),tt(20),tt0(20),CUR_nod(20),CUR_nod_I(20)
      real(kind=kreal) vect(20),vect0(20),d_tt(20) 
      integer(kind=kint) nodLocal(20)
!C
      ALFA = 1.0 - BETA
!
        do k = 1, fstrHEAT%Q_SUF_tot 

          icel    = fstrHEAT%Q_SUF_elem(k)
           jj = hecMESH%global_elem_ID(icel)
          ic_type = hecMESH%elem_type(icel)
          isect   = hecMESH%section_ID(icel)
          IMAT    = hecMESH%section%sect_mat_ID_item(isect)
          isuf    = fstrHEAT%Q_SUF_surf(k)
          iamp    = fstrHEAT%Q_SUF_ampl(k)

          call heat_get_amplitude( fstrHEAT,iamp,CTIME,QQ )
          val     = fstrHEAT%Q_SUF_val (k) * QQ
!C**
          nn = hecmw_get_max_node(ic_type)
!C** 
          iS = hecMESH%elem_node_index(icel-1)
          do j = 1, nn
            nodLOCAL(j) = hecMESH%elem_node_item(iS+j)
            xx(j) = hecMESH%node( 3*nodLOCAL(j)-2 )
            yy(j) = hecMESH%node( 3*nodLOCAL(j)-1 )
            zz(j) = hecMESH%node( 3*nodLOCAL(j)   )
!!!!!    20111014   -----
            tt0 (j) = fstrHEAT%  TEMP0( nodLOCAL(j) )
            tt  (j) = fstrHEAT%  TEMP ( nodLOCAL(j) )
            if ( iterCUR .eq. 1 ) then
              d_tt(j) = fstrHEAT%D_TEMP ( nodLOCAL(j) )
            else if ( iterCUR .gt. 1 ) then
              d_tt(j) = fstrHEAT%D_TEMP ( nodLOCAL(j) )
!             d_tt(j) = fstrHEAT%  TEMP ( nodLOCAL(j) ) - fstrHEAT%TEMP0( nodLOCAL(j) )
            end if

            CUR_nod_I(j) = fstrHEAT%CURE3D_nod_I( nodLOCAL(j)   )
            CUR_nod  (j) = fstrHEAT%CURE3D_nod  ( nodLOCAL(j)   )
!!!!!    20111014   -----
          enddo  
!C**		
          if    ( ic_type.eq.111 ) then
          elseif( ic_type.eq.231 ) then
          elseif( ic_type.eq.232 ) then
          elseif( ic_type.eq.241 ) then
          elseif( ic_type.eq.242 ) then
          elseif( ic_type.eq.341 ) then
          elseif( ic_type.eq.342 ) then
            call cure_DFLUX_342(nn,xx,yy,zz,isuf,val,vect0                            &
                               ,tt0,IMAT,DTIME,icel,jj                                &
                               ,fstrHEAT%RHOtab(IMAT)    ,fstrHEAT%RHOtemp(IMAT,:)    &
                               ,fstrHEAT%RHOfuncA(IMAT,:),fstrHEAT%RHOfuncB(IMAT,:)   &
                               ,fstrHEAT%idx_cure          &
                               ,fstrHEAT%ccc_rl            &
                               ,CUR_nod_I                  &
                               ,IDBG,incr,iterCUR,d_tt)
            call cure_DFLUX_342(nn,xx,yy,zz,isuf,val,vect                             &
                               ,tt,IMAT,DTIME,icel,jj                                 &
                               ,fstrHEAT%RHOtab(IMAT)    ,fstrHEAT%RHOtemp(IMAT,:)    &
                               ,fstrHEAT%RHOfuncA(IMAT,:),fstrHEAT%RHOfuncB(IMAT,:)   &
                               ,fstrHEAT%idx_cure          &
                               ,fstrHEAT%ccc_rl            &
                               ,CUR_nod                    &
                               ,IDBG,incr,iterCUR,d_tt)
          elseif( ic_type.eq.351 ) then
          elseif( ic_type.eq.352 ) then
          elseif( ic_type.eq.361 ) then
          elseif( ic_type.eq.362 ) then
          elseif( ic_type.eq.731 ) then
          elseif( ic_type.eq.741 ) then
          endif
!C 
          do j = 1, nn
!           hecMAT%B( nodLOCAL(j) ) = hecMAT%B( nodLOCAL(j) ) - vect(j) 
            hecMAT%B( nodLOCAL(j) ) = hecMAT%B( nodLOCAL(j) ) - vect(j)*BETA - vect0(j)*ALFA 
          enddo
!C
        enddo
      return
   end subroutine cure_mat_ass_bc_DFLUX 

end module m_heat_mat_ass_bc_DFLUX
