heat_mat_ass_conductivity.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_conductivity
contains
 
  subroutine heat_mat_ass_conductivity( hecMESH,hecMAT,fstrSOLID,fstrHEAT,beta )
    use hecmw
    use m_fstr
    use m_heat_lib
    use m_static_lib
    implicit none
    type(fstr_heat)          :: fstrHEAT
    type(hecmwst_matrix)     :: hecMAT
    type(fstr_solid)         :: fstrSOLID
    type(hecmwst_local_mesh) :: hecMESH
    integer(kind=kint) :: itype, is, iE, ic_type, icel, isect, IMAT, ntab, itab, NDOF
    integer(kind=kint) :: in0, nn, i, in, j, nodLOCAL(20), jsect, ic, ip, inod, jp, jnod, isU, ieU, ik, isL, ieL
    real(kind=kreal)   :: beta, tzero, alpha, temp(1000), funca(1000), funcb(1000), tt(20), t0(20), ss(400)
    real(kind=kreal)   :: asect, thick, gth, ghh, gr1, gr2
    real(kind=kreal) :: lumped(20), stiff(20, 20), ecoord(3,20)
    real(kind=kreal), allocatable :: s(:)
 
    ndof = hecmesh%n_dof
    tzero = hecmesh%zero_temp
    alpha = 1.0 - beta
 
    call hecmw_mat_clear(hecmat)
    call hecmw_mat_clear_b(hecmat)
 
    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)
      ! 761/781 hold the rotational dummy nodes after the translational ones,
      ! and heat analysis uses the translational part only
      if( ic_type == 761 ) ic_type = 731
      if( ic_type == 781 ) ic_type = 741
      if (hecmw_is_etype_link(ic_type)) cycle
      if (hecmw_is_etype_patch(ic_type)) cycle
 
      !$omp parallel default(none), &
        !$omp&  private(icel,isect,IMAT,ntab,itab,temp,funcA,funcB,in0,nn,i,j,nodLOCAL,SS,TT,GTH,GHH,GR1,GR2,ecoord,&
        !$omp&          stiff,in,ASECT,thick,jsect), &
        !$omp&  shared(iS,iE,hecMESH,fstrSOLID,ic_type,hecMAT,fstrHEAT,TZERO)
      !$omp do
      do icel = is, ie
        if( fstrsolid%elements(icel)%elemact_flag == kelact_inactive ) cycle
        isect = hecmesh%section_ID(icel)
        imat = hecmesh%section%sect_mat_ID_item(isect)
 
        if( hecmesh%section%sect_type(isect) .ne. 4 ) then
          ntab = fstrheat%CONDtab(imat)
          do itab = 1, ntab
            temp(itab)  = fstrheat%CONDtemp (imat,itab)
            funca(itab) = fstrheat%CONDfuncA(imat,itab)
            funcb(itab) = fstrheat%CONDfuncB(imat,itab)
          enddo
          funca(ntab+1) = fstrheat%CONDfuncA(imat,ntab+1)
          funcb(ntab+1) = fstrheat%CONDfuncB(imat,ntab+1)
        endif
 
        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)
          tt(i) = fstrheat%TEMP (   nodlocal(i)   )
          !T0(i) = fstrHEAT%TEMP0(   nodLOCAL(i)   )
          do j = 1, 3
            ecoord(j,i) = hecmesh%node(3*(nodlocal(i)-1)+j)
          enddo
        enddo
        do i = 1, nn*nn
          ss(i) = 0.0
        enddo
 
        if(ic_type == 111) then
          in = hecmesh%section%sect_R_index(isect)
          asect = hecmesh%section%sect_R_item(in)
          call heat_conductivity_c1(ic_type, nn, ecoord(1:2,1:nn), tt, imat, asect, stiff, &
            fstrheat%CONDtab(imat), fstrheat%CONDtemp(imat,:), fstrheat%CONDfuncA(imat,:) ,fstrheat%CONDfuncB(imat,:))
 
        elseif(ic_type == 231 .or. ic_type == 232 .or. ic_type == 241 .or. ic_type == 242)then
          in = hecmesh%section%sect_R_index(isect)
          thick = hecmesh%section%sect_R_item(in)
          call heat_conductivity_c2(ic_type, nn, ecoord(1:2,1:nn), tt, imat, thick, stiff, &
            fstrheat%CONDtab(imat), fstrheat%CONDtemp(imat,:), fstrheat%CONDfuncA(imat,:) ,fstrheat%CONDfuncB(imat,:))
 
        elseif(ic_type == 341 .or. ic_type == 342 .or. ic_type == 351 .or. ic_type == 352 .or. &
             & ic_type == 361 .or. ic_type == 362)then
          call heat_conductivity_c3(ic_type, nn, ecoord(1:3,1:nn), tt, imat, stiff, &
            fstrheat%CONDtab(imat), fstrheat%CONDtemp(imat,:), fstrheat%CONDfuncA(imat,:) ,fstrheat%CONDfuncB(imat,:))
 
        elseif (ic_type == 541) then
          jsect = hecmesh%section%sect_R_index(isect-1)+1
          gth = hecmesh%section%sect_R_item(jsect)
          ghh = hecmesh%section%sect_R_item(jsect+1)
          gr1 = hecmesh%section%sect_R_item(jsect+2)
          gr2 = hecmesh%section%sect_R_item(jsect+3)
          call heat_conductivity_541(nn, ecoord, tt, tzero, gth, ghh, gr1, gr2, ss, stiff)
 
        elseif(ic_type == 731)then
          nn = 4
          nodlocal(nn) = hecmesh%elem_node_item(in0+nn-1)
          in = hecmesh%section%sect_R_index(isect)
          thick = hecmesh%section%sect_R_item(in)
          ss = 0.0d0
          call heat_conductivity_shell_731(ic_type, nn, ecoord(1:3,1:nn), tt, imat, thick, ss, stiff, &
            fstrheat%CONDtab(imat), fstrheat%CONDtemp(imat,:), fstrheat%CONDfuncA(imat,:) ,fstrheat%CONDfuncB(imat,:))
 
        elseif(ic_type == 741)then
          in = hecmesh%section%sect_R_index(isect)
          thick = hecmesh%section%sect_R_item(in)
          call heat_conductivity_shell_741(ic_type, nn, ecoord(1:3,1:nn), tt, imat, thick, ss, stiff, &
            fstrheat%CONDtab(imat), fstrheat%CONDtemp(imat,:), fstrheat%CONDfuncA(imat,:) ,fstrheat%CONDfuncB(imat,:))
 
        else
          !$omp critical
          write(*,*)"** error setMASS"
          !$omp end critical
        endif
 
        if(ic_type == 541 .or. ic_type == 731 .or. ic_type == 741)then
          stiff = 0.0d0
          in = 1
          do i = 1, nn
            do j = 1, nn
              stiff(j,i) = ss(in)
              in = in + 1
            enddo
          enddo
        endif
 
        call hecmw_mat_ass_elem(hecmat, nn, nodlocal, stiff)
 
      enddo
      !$omp end do
      !$omp end parallel
    enddo
 
    allocate(s(hecmat%NP))
    s = 0.0d0
 
    call hecmw_matvec(hecmesh, hecmat, fstrheat%TEMP0, s)
 
    hecmat%D  = beta*hecmat%D
    hecmat%AU = beta*hecmat%AU
    hecmat%AL = beta*hecmat%AL
    hecmat%B(1:hecmat%N) = hecmat%B(1:hecmat%N) - alpha*s(1:hecmat%N)
 
    deallocate(s)
 
  end subroutine heat_mat_ass_conductivity
end module m_heat_mat_ass_conductivity