fstr_ass_load.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_fstr_ass_load
  use m_fstr
  use m_static_lib
  use m_fstr_precheck
  use m_fstr_elemact
  use mmechgauss
  use mreadtemp
  use muload
  use m_fstr_spring
  use m_common_struct
  use m_utilities
  
  implicit none
  
contains
  !
  !======================================================================!
 
  subroutine fstr_ass_load(cstep, ctime, hecMESH, hecMAT, fstrSOLID, fstrPARAM)
    !======================================================================!
    integer(kind=kint), intent(in)       :: cstep
    real(kind=kreal), intent(in)         :: ctime  
    type(hecmwst_matrix), intent(inout)  :: hecmat
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    type(fstr_solid), intent(inout)      :: fstrSOLID
    type(fstr_param), intent(inout)      :: fstrPARAM
 
    ! Initialize the global load vector
    fstrsolid%GL(:) = 0.0d0
    fstrsolid%EFORCE(:) = 0.0d0
 
    ! Process concentrated nodal forces (CLOAD)
    call process_concentrated_loads(cstep, ctime, hecmesh, fstrsolid)
    
    ! Process distributed loads (DLOAD) - surface pressure and volume force
    call process_distributed_loads(cstep, ctime, hecmesh, fstrsolid)
    
    ! Process user-defined loads
    call process_user_loads(cstep, fstrsolid)
    
    ! Update global load vector
    call hecmw_update_r(hecmesh, fstrsolid%GL, hecmesh%n_node, hecmesh%n_dof)
    
    ! Update right-hand side vector
    call hecmw_mat_clear_b(hecmat)
    call update_rhs_vector(hecmesh, hecmat, fstrsolid)
    
    ! Process thermal loads (TLOAD)
    call process_thermal_loads(cstep, ctime, hecmesh, hecmat, fstrsolid, fstrparam)
    
    ! Process spring forces
    call fstr_update_ndforce_spring(cstep, hecmesh, fstrsolid, hecmat%B)
    
  end subroutine fstr_ass_load
 
  !======================================================================!
  !======================================================================!
  subroutine process_concentrated_loads(cstep, ctime, hecMESH, fstrSOLID)
    integer(kind=kint), intent(in)       :: cstep
    real(kind=kreal), intent(in)         :: ctime
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    type(fstr_solid), intent(inout)      :: fstrSOLID
    
    integer(kind=kint) :: n_rot, rid, n_nodes, idof, ndof
    integer(kind=kint) :: ig0, ig, ityp, iS0, iE0, ik, in, grpid, jj_n_amp
    real(kind=kreal)   :: aval, fval, tval
    real(kind=kreal)   :: normal(3), ccoord(3), cdisp(3), cdiff(3)
    real(kind=kreal)   :: vect(60)
    type(trotinfo)     :: rinfo
    
    ndof = hecmesh%n_dof
 
    ! Initialize rotation information for torque loads
    n_rot = fstrsolid%CLOAD_ngrp_rot
    if (n_rot > 0) call fstr_rotinfo_init(n_rot, rinfo)
 
    ! Process all concentrated loads
    do ig0 = 1, fstrsolid%CLOAD_ngrp_tot
      grpid = fstrsolid%CLOAD_ngrp_GRPID(ig0)
      if (.not. fstr_isloadactive(fstrsolid, grpid, cstep)) cycle
      ig = fstrsolid%CLOAD_ngrp_ID(ig0)
      ityp = fstrsolid%CLOAD_ngrp_DOF(ig0)
      fval = fstrsolid%CLOAD_ngrp_val(ig0)
      jj_n_amp = fstrsolid%CLOAD_ngrp_amp(ig0)
      aval = fval
      if (jj_n_amp <= 0) then  ! Amplitude not defined
        aval = fstrsolid%FACTOR(2) * fval
      else
        call fstr_get_amplitude(hecmesh, fstrsolid, cstep, jj_n_amp, ctime, aval)
      endif
 
      ! A load continuing from a previous step is held at its nominal value,
      ! BUT only when no amplitude is given. An explicit AMP= prescribes the
      ! load's time history, so it is honored in every active step (consistent
      ! with the dynamic load path, which has no such override).
      if (jj_n_amp <= 0 .and. fstr_isloadactive(fstrsolid, grpid, cstep-1)) aval = fval
      is0 = hecmesh%node_group%grp_index(ig-1) + 1
      ie0 = hecmesh%node_group%grp_index(ig)
 
      if( fstrsolid%CLOAD_ngrp_rotID(ig0) > 0 ) then ! setup torque load information
        rid = fstrsolid%CLOAD_ngrp_rotID(ig0)
        if (.not. rinfo%conds(rid)%active) then
          rinfo%conds(rid)%active = .true.
          rinfo%conds(rid)%center_ngrp_id = fstrsolid%CLOAD_ngrp_centerID(ig0)
          rinfo%conds(rid)%torque_ngrp_id = ig
        endif
        if (ityp > ndof) ityp = ityp - ndof
        rinfo%conds(rid)%vec(ityp) = aval
        cycle
      endif
 
      do ik = is0, ie0
        in = hecmesh%node_group%grp_item(ik)
        fstrsolid%GL(ndof*(in-1)+ityp) = fstrsolid%GL(ndof*(in-1)+ityp) + aval
      enddo
    enddo
 
    !Add torque load to fstrSOLID%GL
    do rid = 1, n_rot
      if (.not. rinfo%conds(rid)%active) cycle
      ! Get number of slave nodes
      n_nodes = hecmw_ngrp_get_number(hecmesh, rinfo%conds(rid)%torque_ngrp_id)
 
      ! Get center node
      ig = rinfo%conds(rid)%center_ngrp_id
      do idof = 1, ndof
        ccoord(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, hecmesh%node)
        cdisp(idof) = hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, fstrsolid%unode)
        cdisp(idof) = cdisp(idof) + hecmw_ngrp_get_totalvalue(hecmesh, ig, ndof, idof, fstrsolid%dunode)
      enddo
      ccoord(1:ndof) = ccoord(1:ndof) + cdisp(1:ndof)
 
      tval = dsqrt(dot_product(rinfo%conds(rid)%vec(1:ndof), rinfo%conds(rid)%vec(1:ndof)))
      if (tval < 1.d-16) then
        write(*,*) '###ERROR### : norm of torque vector must be > 0.0'
        call hecmw_abort(hecmw_comm_get_comm())
      endif
      normal(1:ndof) = rinfo%conds(rid)%vec(1:ndof) / tval
      tval = tval / dble(n_nodes)
 
      ig = rinfo%conds(rid)%torque_ngrp_id
      is0 = hecmesh%node_group%grp_index(ig-1) + 1
      ie0 = hecmesh%node_group%grp_index(ig)
      do ik = is0, ie0
        in = hecmesh%node_group%grp_item(ik)
        cdiff(1:ndof) = hecmesh%node(ndof*(in-1)+1:ndof*in) + fstrsolid%unode(ndof*(in-1)+1:ndof*in) &
          & + fstrsolid%dunode(ndof*(in-1)+1:ndof*in) - ccoord(1:ndof)
        call cross_product(normal,cdiff,vect(1:ndof))
        fval = dot_product(vect(1:ndof), vect(1:ndof))
        if (fval < 1.d-16) then
          write(*,*) '###ERROR### : torque node is at the same position as that of center node in rotational surface.'
          call hecmw_abort(hecmw_comm_get_comm())
        endif
        vect(1:ndof) = (tval/fval) * vect(1:ndof)
        fstrsolid%GL(ndof*(in-1)+1:ndof*in) = fstrsolid%GL(ndof*(in-1)+1:ndof*in) + vect(1:ndof)
      enddo
    enddo
    if (n_rot > 0) call fstr_rotinfo_finalize(rinfo)
    
  end subroutine process_concentrated_loads
    !
    ! -------------------------------------------------------------------
    !  DLOAD
    ! -------------------------------------------------------------------
  subroutine process_distributed_loads(cstep, ctime, hecMESH, fstrSOLID)
    integer(kind=kint), intent(in)       :: cstep
    real(kind=kreal), intent(in)         :: ctime
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    type(fstr_solid), intent(inout)      :: fstrSOLID
    
    integer(kind=kint) :: ndof, ig0, ig, ltype, iS0, iE0, ik, icel, ic_type, nn, is
    integer(kind=kint) :: isect, id, iset, ihead, nsize, grpid, i, j, jj_n_amp
    integer(kind=kint) :: iwk(60), nodLocal(20)
    real(kind=kreal)   :: xx(20), yy(20), zz(20), vect(60), params(0:6)
    real(kind=kreal)   :: factor, aval, rho, thick, pa1
    logical            :: fg_surf
    type(tmaterial), pointer :: material
    
    ndof = hecmesh%n_dof
    
    do ig0 = 1, fstrsolid%DLOAD_ngrp_tot
      grpid = fstrsolid%DLOAD_ngrp_GRPID(ig0)
      if (.not. fstr_isloadactive(fstrsolid, grpid, cstep)) cycle
      jj_n_amp = fstrsolid%DLOAD_ngrp_amp(ig0)
      ig = fstrsolid%DLOAD_ngrp_ID(ig0)
      ltype = fstrsolid%DLOAD_ngrp_LID(ig0)
      do i = 0, 6
        params(i) = fstrsolid%DLOAD_ngrp_params(i, ig0)
      enddo
      if (jj_n_amp <= 0) then  ! Amplitude not defined
        factor = fstrsolid%factor(2)
      else
        ! params(0) is the load magnitude (pressure / body force value,
        ! gravity acceleration, or angular velocity). Passing it to the
        ! evaluator lets VALUE=ABSOLUTE replace the magnitude with the
        ! amplitude value, while VALUE=RELATIVE reduces the factor to a(t).
        aval = params(0)
        call fstr_get_amplitude(hecmesh, fstrsolid, cstep, jj_n_amp, ctime, aval)
        if (dabs(params(0)) > 1.d-30) then
          factor = aval / params(0)
        else
          factor = 1.0d0
        endif
      endif
 
      ! Continuing load held at nominal only when no amplitude is given;
      ! an explicit AMP= is honored in every active step (see CLOAD above).
      if (jj_n_amp <= 0 .and. fstr_isloadactive(fstrsolid, grpid, cstep-1)) factor = 1.0d0
      ! ----- START & END
      fg_surf = (ltype == 100)
      if( fg_surf ) then                  ! surface group
        is0 = hecmesh%surf_group%grp_index(ig-1) + 1
        ie0 = hecmesh%surf_group%grp_index(ig)
      else                                ! element group
        is0 = hecmesh%elem_group%grp_index(ig-1) + 1
        ie0 = hecmesh%elem_group%grp_index(ig)
      endif
      do ik = is0, ie0
        if( fg_surf ) then                ! surface group
          ltype = hecmesh%surf_group%grp_item(2*ik) * 10
          icel = hecmesh%surf_group%grp_item(2*ik-1)
          ic_type = hecmesh%elem_type(icel)
        else                              ! element group
          icel = hecmesh%elem_group%grp_item(ik)
          ic_type = hecmesh%elem_type(icel)
        endif
 
        !ELEMENT ACTIVATION
        if( fstrsolid%elements(icel)%elemact_flag == kelact_inactive ) cycle
 
        if (hecmw_is_etype_link(ic_type)) cycle
        if (hecmw_is_etype_patch(ic_type)) cycle
        ! if( ic_type==3422 ) ic_type=342
        nn = hecmw_get_max_node(ic_type)
        ! ----- node ID
        is = hecmesh%elem_node_index(icel-1)
        if (fstrsolid%DLOAD_follow == 0) then
          do j = 1, nn
            nodlocal(j) = hecmesh%elem_node_item (is+j)
            ! ----- 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)   )
            ! ----- create iwk array ***
            do i = 1, ndof
              iwk( ndof*(j-1)+i ) = ndof*( nodlocal(j)-1 )+i
            enddo
          enddo
        else
          do j = 1, nn
            nodlocal(j) = hecmesh%elem_node_item (is+j)
            ! ----- nodal coordinate
            if (ndof==2) then
              xx(j) = hecmesh%node( 3*nodlocal(j)-2 )+fstrsolid%unode( 2*nodlocal(j)-1 )+fstrsolid%dunode( 2*nodlocal(j)-1 )
              yy(j) = hecmesh%node( 3*nodlocal(j)-1 )+fstrsolid%unode( 2*nodlocal(j)   )+fstrsolid%dunode( 2*nodlocal(j)   )
            else if (ndof==3) then
              xx(j) = hecmesh%node( 3*nodlocal(j)-2 )+fstrsolid%unode( 3*nodlocal(j)-2 )+fstrsolid%dunode( 3*nodlocal(j)-2 )
              yy(j) = hecmesh%node( 3*nodlocal(j)-1 )+fstrsolid%unode( 3*nodlocal(j)-1 )+fstrsolid%dunode( 3*nodlocal(j)-1 )
              zz(j) = hecmesh%node( 3*nodlocal(j)   )+fstrsolid%unode( 3*nodlocal(j)   )+fstrsolid%dunode( 3*nodlocal(j)   )
            else if (ndof==6) then
              xx(j) = hecmesh%node( 3*nodlocal(j)-2 )+fstrsolid%unode( 6*nodlocal(j)-5 )+fstrsolid%dunode( 6*nodlocal(j)-5 )
              yy(j) = hecmesh%node( 3*nodlocal(j)-1 )+fstrsolid%unode( 6*nodlocal(j)-4 )+fstrsolid%dunode( 6*nodlocal(j)-4 )
              zz(j) = hecmesh%node( 3*nodlocal(j)   )+fstrsolid%unode( 6*nodlocal(j)-3 )+fstrsolid%dunode( 6*nodlocal(j)-3 )
            endif
            ! ----- create iwk array ***
            do i = 1, ndof
              iwk( ndof*(j-1)+i ) = ndof*( nodlocal(j)-1 )+i
            enddo
          enddo
        endif
        ! ----- section  ID
        isect = hecmesh%section_ID(icel)
        ! ----- Get Properties
        material => fstrsolid%elements(icel)%gausses(1)%pMaterial
        rho = material%variables(m_density)
        call fstr_get_thickness(hecmesh, isect, thick)
        ! ----- Section Data
        if (ndof == 2) then
          id = hecmesh%section%sect_opt(isect)
          if (id == 0) then
            iset = 1
          else if (id == 1) then
            iset = 0
          else if (id == 2) then
            iset = 2
          endif
          pa1 = 1.d0
        endif
        ! ----- Create local stiffness
        if (ic_type==301)then
          ihead = hecmesh%section%sect_R_index(isect-1)
          call dl_c1(nn,xx(1:nn),yy(1:nn),zz(1:nn),rho,thick,ltype,params,vect(1:nn*ndof),nsize)
 
        elseif( ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 .or. ic_type == 2322 ) then
          call dl_c2(ic_type,nn,xx(1:nn),yy(1:nn),rho,pa1,ltype,params,vect(1:nn*ndof),nsize,iset)
 
        else if ( ic_type == 341 .or. ic_type == 351 .or. ic_type == 361 .or.   &
            ic_type == 342 .or. ic_type == 352 .or. ic_type == 362 ) then
          call dl_c3(ic_type,nn,xx(1:nn),yy(1:nn),zz(1:nn),rho,ltype,params,vect(1:nn*ndof),nsize)
 
        else if ( ic_type == 641 ) then
          ihead = hecmesh%section%sect_R_index(isect-1)
          call dl_beam_641(ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), rho, ltype, params, &
            hecmesh%section%sect_R_item(ihead+1:), vect(1:nn*ndof), nsize)
 
        else if( ( ic_type == 741 ) .or. ( ic_type == 743 ) .or. ( ic_type == 731 ) ) then
          call dl_shell(ic_type, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, fstrsolid%elements(icel)%gausses)
 
        else if( ( ic_type==761 ) .or. ( ic_type==781 ) ) then
          call dl_shell_33(ic_type, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, &
            fstrsolid%elements(icel)%gausses)
 
        else
          nsize = 0
          write(*,*)"### WARNING: DLOAD",ic_type
 
        endif
        ! ----- Add vector
        do j = 1, nsize
          fstrsolid%GL(iwk(j)) = fstrsolid%GL(iwk(j)) + factor * vect(j)
        enddo
      enddo
    enddo
    
  end subroutine process_distributed_loads
 
    ! -----Uload
  subroutine process_user_loads(cstep, fstrSOLID)
    integer(kind=kint), intent(in)  :: cstep
    type(fstr_solid), intent(inout) :: fstrSOLID
    
    real(kind=kreal) :: factor
    
    factor = fstrsolid%factor(2)
    call uloading(cstep, factor, fstrsolid%GL)
    
  end subroutine process_user_loads
 
  !======================================================================!
  !======================================================================!
  subroutine update_rhs_vector(hecMESH, hecMAT, fstrSOLID)
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    type(hecmwst_matrix), intent(inout)  :: hecMAT
    type(fstr_solid), intent(inout)      :: fstrSOLID
    
    integer(kind=kint) :: i
    
    do i = 1, hecmesh%n_node * hecmesh%n_dof
      hecmat%B(i) = fstrsolid%GL(i) - fstrsolid%QFORCE(i)
    enddo
    
    do i = 1, hecmat%NDOF * hecmat%NP
      !thermal load is not considered
      fstrsolid%EFORCE(i) = fstrsolid%GL(i)
    enddo
    
  end subroutine update_rhs_vector
 
    ! -------------------------------------------------------------------
    !  TLOAD : THERMAL LOAD USING TEMPERATURE
    ! -------------------------------------------------------------------
  subroutine process_thermal_loads(cstep, ctime, hecMESH, hecMAT, fstrSOLID, fstrPARAM)
    integer(kind=kint), intent(in)       :: cstep
    real(kind=kreal), intent(in)         :: ctime
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    type(hecmwst_matrix), intent(inout)  :: hecMAT
    type(fstr_solid), intent(inout)      :: fstrSOLID
    type(fstr_param), intent(in)         :: fstrPARAM
    
    integer(kind=kint) :: ndof, ig0, ig, iS0, iE0, ik, in, grpid
    integer(kind=kint) :: itype, is, iE, icel, ic_type, nn, isect, cdsys_ID, id, iset
    integer(kind=kint) :: i, j, ihead, nodLocal(20), iwk(60)
    real(kind=kreal)   :: factor, fval, pa1
    real(kind=kreal)   :: xx(20), yy(20), zz(20), tt(20), tt0(20), coords(3,3), vect(60)
    real(kind=kreal)   :: local_coords(3,3)  ! Local copy for coordinate transformation
    
    ndof = hecmesh%n_dof
    
    if (fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0) then
      do ig0 = 1, fstrsolid%TEMP_ngrp_tot
        grpid = fstrsolid%TEMP_ngrp_GRPID(ig0)
        if (.not. fstr_isloadactive(fstrsolid, grpid, cstep)) cycle
        factor = fstrsolid%factor(2)
        if (fstr_isloadactive(fstrsolid, grpid, cstep-1)) factor = 1.0d0
        ig = fstrsolid%TEMP_ngrp_ID(ig0)
        fval = fstrsolid%TEMP_ngrp_val(ig0)
        is0 = hecmesh%node_group%grp_index(ig-1) + 1
        ie0 = hecmesh%node_group%grp_index(ig)
        do ik = is0, ie0
          in = hecmesh%node_group%grp_item(ik)
          pa1 = fstrsolid%temp_bak(in)
          fstrsolid%temperature(in) = pa1 + (fval - pa1) * factor
        enddo
      enddo
      
      if (fstrsolid%TEMP_irres > 0) then
        call read_temperature_result(hecmesh, fstrsolid%TEMP_irres, fstrsolid%TEMP_tstep, &
          &  fstrsolid%TEMP_rtype, fstrsolid%TEMP_interval, fstrsolid%TEMP_factor, ctime, &
          &  fstrsolid%temperature, fstrsolid%temp_bak)
      endif
    endif
 
    ! ----- elemact element
    if( fstrsolid%elemact%ELEMACT_egrp_tot > 0 ) &
      &  call fstr_update_elemact_solid( hecmesh, fstrsolid, cstep, ctime )
 
    if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 ) then
      ! ----- element TYPE loop.
      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
        ! ----- Set number of nodes
        nn = hecmw_get_max_node(ic_type)
        
        ! ----- element loop
        do icel = is, ie
 
          !ELEMENT ACTIVATION
          if( fstrsolid%elements(icel)%elemact_flag == kelact_inactive ) cycle
 
          ! ----- node ID
          is = hecmesh%elem_node_index(icel-1)
          do j = 1, nn
            nodlocal(j) = hecmesh%elem_node_item(is+j)
            ! ----- nodal coordinate
            if (ndof == 2) then
              xx(j) = hecmesh%node(3*nodlocal(j)-2) + fstrsolid%unode(ndof*nodlocal(j)-1)
              yy(j) = hecmesh%node(3*nodlocal(j)-1) + fstrsolid%unode(ndof*nodlocal(j))
            else if (ndof == 3) then
              xx(j) = hecmesh%node(3*nodlocal(j)-2) + fstrsolid%unode(ndof*nodlocal(j)-2)
              yy(j) = hecmesh%node(3*nodlocal(j)-1) + fstrsolid%unode(ndof*nodlocal(j)-1)
              zz(j) = hecmesh%node(3*nodlocal(j)) + fstrsolid%unode(ndof*nodlocal(j))
            endif
            tt0(j) = fstrsolid%last_temp(nodlocal(j))
            tt(j) = fstrsolid%temperature(nodlocal(j))
            ! ----- create iwk array ***
            do i = 1, ndof
              iwk(ndof*(j-1)+i) = ndof*(nodlocal(j)-1)+i
            enddo
          enddo
          
          ! ----- section  Data
          isect = hecmesh%section_ID(icel)
          cdsys_id = hecmesh%section%sect_orien_ID(isect)
          call get_coordsys(cdsys_id, hecmesh, fstrsolid, coords, icel)
 
          if( ndof == 2 ) then
            id=hecmesh%section%sect_opt(isect)
            if( id==0 ) then
              iset = 1
            else if (id == 1) then
              iset = 0
            else if (id == 2) then
              iset = 2
            endif
            pa1 = 1.0d0
          endif
          
          if (ic_type == 641) then
            isect = hecmesh%section_ID(icel)
            ihead = hecmesh%section%sect_R_index(isect-1)
            
            call tload_beam_641( nn, ndof, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn),    &
              fstrsolid%elements(icel)%gausses, hecmesh%section%sect_R_item(ihead+1:), &
              fstrparam%ref_temp, &
              vect(1:nn*ndof))
            
            do j = 1, ndof*nn
              hecmat%B(iwk(j)) = hecmat%B(iwk(j)) + vect(j)
            enddo
            cycle
          endif
          
          ! Local copy of coordinate data
          local_coords = coords
          
          ! Calculate thermal load based on element type
          call calculate_thermal_load(ic_type, nn, xx, yy, zz, tt, tt0, isect, ndof, &
                                     hecmesh, fstrsolid, fstrparam, icel, vect, cdsys_id, local_coords, &
                                     iset, pa1, iwk, hecmat%B)
        enddo
      enddo
    endif
    
  end subroutine process_thermal_loads
 
  !======================================================================!
  !======================================================================!
  subroutine calculate_thermal_load(ic_type, nn, xx, yy, zz, tt, tt0, isect, ndof, &
                                   hecMESH, fstrSOLID, fstrPARAM, icel, vect, cdsys_ID, coords, &
                                   iset, pa1, iwk, B)
    integer(kind=kint), intent(in)       :: ic_type, nn, isect, ndof, cdsys_ID, iset
    real(kind=kreal), intent(in)         :: xx(*), yy(*), zz(*), tt(*), tt0(*), pa1
    real(kind=kreal), intent(inout)      :: coords(3,3)  ! Changed INTENT from IN to INOUT
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    type(fstr_solid), intent(in)         :: fstrSOLID
    type(fstr_param), intent(in)         :: fstrPARAM
    integer(kind=kint), intent(in)       :: icel
    real(kind=kreal), intent(out)        :: vect(*)
    integer(kind=kint), intent(in)       :: iwk(*)
    real(kind=kreal), intent(inout)      :: b(*)
    
    integer(kind=kint) :: j, myrank
    
    myrank = 0
    
    ! 2D elements
    if (ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232) then
      call tload_c2(ic_type, nn, xx(1:nn), yy(1:nn), tt(1:nn), tt0(1:nn), &
                   fstrsolid%elements(icel)%gausses, pa1, iset, vect(1:nn*2))
    
    else if (ic_type == 361) then
      if (fstrsolid%sections(isect)%elemopt361 == kel361fi) then
              call tload_c3                                                          &
                ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), fstrparam%ref_temp, &
                     fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords)
      else if (fstrsolid%sections(isect)%elemopt361 == kel361bbar) then
              call tload_c3d8bbar                                                          &
                ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), fstrparam%ref_temp, &
                           fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords)
      else if (fstrsolid%sections(isect)%elemopt361 == kel361ic) then
              call tload_c3d8ic                                                            &
                ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), fstrparam%ref_temp, &
                         fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords)
      else if (fstrsolid%sections(isect)%elemopt361 == kel361fbar) then
              call tload_c3d8fbar                                                            &
                ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), fstrparam%ref_temp, &
                           fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords)
      else if (fstrsolid%sections(isect)%elemopt361 == kel361up) then
              ! UP element : reuse the standard C3 thermal load. Note the thermal
              ! strain is volumetric and TLOAD_C3 uses the full elastic matrix, so
              ! this is not strictly consistent with the u-p deviatoric/pressure
              ! split (see report); provided for basic thermal loading support.
              call tload_c3                                                          &
                ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), fstrparam%ref_temp, &
                fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords )
      endif
    
    else if (ic_type == 341 .or. ic_type == 351 .or. &
             ic_type == 342 .or. ic_type == 352 .or. ic_type == 362) then
            call tload_c3                                                                &
              ( ic_type, nn, xx(1:nn), yy(1:nn), zz(1:nn), tt(1:nn), tt0(1:nn), fstrparam%ref_temp, &
                   fstrsolid%elements(icel)%gausses, vect(1:nn*ndof), cdsys_id, coords)
    
    else if (ic_type == 741 .or. ic_type == 743 .or. ic_type == 731 .or. &
             ic_type == 761 .or. ic_type == 781) then
      if (myrank == 0) then
        write(imsg,*) '*------------------------', &
                     '-------------------*'
        write(imsg,*) ' Thermal loading option for shell elements', &
                     'not yet available.'
        write(imsg,*) '*------------------------', &
                     '-------------------*'
        call hecmw_abort(hecmw_comm_get_comm())
      endif
    endif
    
          ! ----- Add vector
    do j = 1, ndof*nn
      b(iwk(j)) = b(iwk(j)) + vect(j)
    enddo
    
  end subroutine calculate_thermal_load
 
end module m_fstr_ass_load