heat_mat_ass_bc_FILM.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
module m_heat_mat_ass_bc_film
contains
  !C
  !C***
  !C*** MAT_ASS_FILM
  !C***
  !C
  subroutine heat_mat_ass_bc_film( hecMESH, hecMAT, fstrSOLID, fstrHEAT, CTIME )
 
    use m_fstr
    use m_heat_get_amplitude
    use m_heat_lib_film
 
    implicit none
    integer(kind=kint) :: k, icel, isuf, iam1, iam2, ic_type, isect, nn, is, j, mm, m, ic, ip
    integer(kind=kint) :: inod, jp, jnod, isU, ieU, ik, isL, ieL
    real(kind=kreal)   :: ctime, qq, hh, sink, thick
    type(fstr_solid)         :: fstrsolid
    type(fstr_heat)          :: fstrHEAT
    type(hecmwst_matrix)     :: hecMAT
    type(hecmwst_local_mesh) :: hecMESH
 
    real(kind=kreal)   :: xx(20), yy(20), zz(20)
    real(kind=kreal)   :: term1(64), term2(20), stiff(8,8)
    integer(kind=kint) :: nodLocal(20), nsuf(8), nodSurf(8)
 
!C
    !$omp parallel default(none), &
      !$omp&  private(k,icel,isuf,iam1,iam2,QQ,HH,SINK,ic_type,isect,nn,is,j,nodLocal, &
      !$omp&  xx,yy,zz,thick,mm,term1,term2,stiff,nsuf,nodSurf,ip,inod,jnod,ic,isU,ieU,ik,jp,isL,ieL), &
      !$omp&  shared(fstrSOLID,fstrHEAT,CTIME,hecMAT,hecMESH)
    !$omp do
    do k = 1, fstrheat%H_SUF_tot
      icel    = fstrheat%H_SUF_elem(k)
      if( fstrsolid%elements(icel)%elemact_flag == kelact_inactive ) cycle
 
      isuf    = fstrheat%H_SUF_surf(k)
      iam1    = fstrheat%H_SUF_ampl(k,1)
      iam2    = fstrheat%H_SUF_ampl(k,2)
      call heat_get_amplitude ( fstrheat, iam1, ctime, qq )
      hh      = fstrheat%H_SUF_val (k,1) * qq
      if( dabs(hh) < 1.d-16 ) cycle
      call heat_get_amplitude ( fstrheat, iam2, ctime, qq )
      sink    = fstrheat%H_SUF_val (k,2) * qq
      ic_type = hecmesh%elem_type(icel)
      if( ic_type == 761 ) ic_type = 731
      if( ic_type == 781 ) ic_type = 741
      isect   = hecmesh%section_ID(icel)
      !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)   )
      enddo
      !C**
      if    ( ic_type.eq.231 ) then
        is = hecmesh%section%sect_R_index(isect)
        thick = hecmesh%section%sect_R_item(is)
        mm = 2
        call heat_film_231(nn,xx,yy,thick,isuf,hh,sink,mm,term1,term2,nsuf)
      elseif( ic_type.eq.232 ) then
        is = hecmesh%section%sect_R_index(isect)
        thick = hecmesh%section%sect_R_item(is)
        mm = 3
        call heat_film_232(nn,xx,yy,thick,isuf,hh,sink,mm,term1,term2,nsuf)
      elseif( ic_type.eq.241 ) then
        is = hecmesh%section%sect_R_index(isect)
        thick = hecmesh%section%sect_R_item(is)
        mm = 2
        call heat_film_241(nn,xx,yy,thick,isuf,hh,sink,mm,term1,term2,nsuf)
      elseif( ic_type.eq.242 ) then
        is = hecmesh%section%sect_R_index(isect)
        thick = hecmesh%section%sect_R_item(is)
        mm = 3
        call heat_film_242(nn,xx,yy,thick,isuf,hh,sink,mm,term1,term2,nsuf)
      elseif( ic_type.eq.341 ) then
        mm = 3
        call heat_film_341(nn,xx,yy,zz,isuf,hh,sink,mm,term1,term2,nsuf)
      elseif( ic_type.eq.342 ) then
        mm = 6
        call heat_film_342(nn,xx,yy,zz,isuf,hh,sink,mm,term1,term2,nsuf)
      elseif( ic_type.eq.351 ) then
        mm = 4
        if( isuf.eq.1 .or. isuf.eq.2 ) mm = 3
        call heat_film_351(nn,xx,yy,zz,isuf,hh,sink,mm,term1,term2,nsuf)
      elseif( ic_type.eq.352 ) then
        mm = 8
        if( isuf.eq.1 .or. isuf.eq.2 ) mm = 6
        call heat_film_352(nn,xx,yy,zz,isuf,hh,sink,mm,term1,term2,nsuf)
      elseif( ic_type.eq.361 ) then
        mm = 4
        call heat_film_361(nn,xx,yy,zz,isuf,hh,sink,mm,term1,term2,nsuf)
      elseif( ic_type.eq.362 ) then
        mm = 8
        call heat_film_362(nn,xx,yy,zz,isuf,hh,sink,mm,term1,term2,nsuf)
      elseif( ic_type.eq.731 ) then
        mm = 3
        call heat_film_731(nn,xx,yy,zz,hh,sink,term1,term2)
        do m = 1, mm
          nsuf(m) = m
        enddo
 
      elseif( ic_type.eq.741 ) then
        mm = 4
        call heat_film_741(nn,xx,yy,zz,hh,sink,term1,term2)
        do m = 1, mm
          nsuf(m) = m
        enddo
 
      endif
      !C
 
      do ip = 1, mm
        nodsurf(ip) = nodlocal(nsuf(ip))
      end do
 
      ic = 0
      stiff = 0.d0
      do ip = 1, mm
        do jp = 1, mm
          ic = ic + 1
          stiff(ip,jp) = -term1(ic)
        enddo
      enddo
 
      call hecmw_mat_ass_elem(hecmat, mm, nodsurf, stiff)
 
      do ip = 1, mm
        !$omp atomic
        hecmat%B(nodsurf(ip)) = hecmat%B(nodsurf(ip)) - term2(ip)
      end do
 
      !C
      !C
    enddo
    !$omp end do
    !$omp end parallel
 
  end subroutine  heat_mat_ass_bc_film
end module m_heat_mat_ass_bc_film