heat_mat_ass_capacity.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_capacity
contains
 
  subroutine heat_mat_ass_capacity(hecMESH, hecMAT, fstrSOLID, fstrHEAT, delta_time)
    use m_fstr
    use m_heat_lib_capacity
    implicit none
    type(fstr_heat)          :: fstrHEAT
    type(hecmwst_matrix)     :: hecMAT
    type(fstr_solid)         :: fstrSOLID
    type(hecmwst_local_mesh) :: hecMESH
    integer(kind=kint) :: i, in, j, nodLOCAL(20), ip, inod
    integer(kind=kint) :: itype, iS, iE, ic_type, icel, isect, IMAT, in0, nn
    real(kind=kreal) :: temp(20), lumped(120), mass(20*6, 20*6), ecoord(3,20)
    real(kind=kreal) :: delta_time, surf, thick, alfa, beta
 
    beta = fstrheat%beta
 
    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 (hecmw_is_etype_link(ic_type)) cycle
      if (hecmw_is_etype_patch(ic_type)) cycle
      if(ic_type == 3414) cycle
 
      !$omp parallel default(none), &
        !$omp&  private(icel,isect,IMAT,nn,temp,in0,i,j,nodLOCAL,ecoord,in,thick,surf,inod,lumped,mass), &
        !$omp&  shared(iS,iE,hecMESH,fstrSOLID,ic_type,hecMAT,fstrHEAT,delta_time)
      !$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)
        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)
          do j = 1, 3
            ecoord(j,i) = hecmesh%node(3*(nodlocal(i)-1)+j)
          enddo
          temp(i) = fstrheat%TEMP0(nodlocal(i))
        enddo
 
        lumped = 0.0d0
        if(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_capacity_c2(ic_type, nn, ecoord(1:2,1:nn), temp, imat, thick, lumped, mass, &
            fstrheat%CPtab(imat), fstrheat%CPtemp(imat,:), fstrheat%CPfuncA(imat,:) ,fstrheat%CPfuncB(imat,:), &
            fstrheat%RHOtab(imat), fstrheat%RHOtemp(imat,:), fstrheat%RHOfuncA(imat,:), fstrheat%RHOfuncB(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_capacity_c3(ic_type, nn, ecoord(1:3,1:nn), temp, imat, lumped, mass, &
            fstrheat%CPtab(imat), fstrheat%CPtemp(imat,:), fstrheat%CPfuncA(imat,:) ,fstrheat%CPfuncB(imat,:), &
            fstrheat%RHOtab(imat), fstrheat%RHOtemp(imat,:), fstrheat%RHOfuncA(imat,:), fstrheat%RHOfuncB(imat,:))
 
        elseif (ic_type == 541) then
          ! skip warning
 
        elseif(ic_type == 731)then
          in = hecmesh%section%sect_R_index(isect)
          thick = hecmesh%section%sect_R_item(in)
          call heat_capacity_shell_731(ic_type, nn, ecoord(1:3,1:nn), temp, imat, thick, lumped, mass, &
            fstrheat%CPtab(imat), fstrheat%CPtemp(imat,:), fstrheat%CPfuncA(imat,:) ,fstrheat%CPfuncB(imat,:), &
            fstrheat%RHOtab(imat), fstrheat%RHOtemp(imat,:), fstrheat%RHOfuncA(imat,:), fstrheat%RHOfuncB(imat,:))
 
        elseif(ic_type == 741)then
          in = hecmesh%section%sect_R_index(isect)
          thick = hecmesh%section%sect_R_item(in)
          call heat_capacity_shell_741(ic_type, nn, ecoord(1:3,1:nn), temp, imat, thick, lumped, mass, &
            fstrheat%CPtab(imat), fstrheat%CPtemp(imat,:), fstrheat%CPfuncA(imat,:) ,fstrheat%CPfuncB(imat,:), &
            fstrheat%RHOtab(imat), fstrheat%RHOtemp(imat,:), fstrheat%RHOfuncA(imat,:), fstrheat%RHOfuncB(imat,:))
 
        elseif(ic_type == 111 .or. ic_type == 301 .or. ic_type == 611 .or. ic_type == 641)then
          in = hecmesh%section%sect_R_index(isect)
          surf = hecmesh%section%sect_R_item(in)
          call heat_capacity_c1(ic_type, nn, ecoord(1:3,1:nn), temp, imat, surf, lumped, mass, &
            fstrheat%CPtab(imat), fstrheat%CPtemp(imat,:), fstrheat%CPfuncA(imat,:) ,fstrheat%CPfuncB(imat,:), &
            fstrheat%RHOtab(imat), fstrheat%RHOtemp(imat,:), fstrheat%RHOfuncA(imat,:), fstrheat%RHOfuncB(imat,:))
 
        else
          write(*,*)"** error setMASS"
        endif
 
        do ip = 1, nn
          inod = nodlocal(ip)
          !$omp atomic
          hecmat%D(inod) = hecmat%D(inod) + lumped(ip) / delta_time
          !$omp atomic
          hecmat%B(inod) = hecmat%B(inod) + lumped(ip)*temp(ip) / delta_time
        enddo
      enddo
      !$omp end do
      !$omp end parallel
    enddo
  end subroutine heat_mat_ass_capacity
end module m_heat_mat_ass_capacity