fstr_frequency_analysis.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 fstr_frequency_analysis
 
  use m_fstr
  use m_fstr_creatematrix_and_dampingforce
  use m_fstr_addbc
  use m_fstr_eig_setmass
  use m_fstr_nodalkinematics
  use m_fstr_update
  use m_fstr_residual, only: fstr_update_reaction_spc
  use m_fstr_nodalstress
  use m_make_result
  use m_hecmw2fstr_mesh_conv
  use hecmw_setup_util, only: node_global_to_local
  use hecmw_ebc_defer
 
  implicit none
 
contains
 
  subroutine fstr_solve_frequency_analysis(hecMESH, hecMAT, fstrSOLID, fstrEIG,  &
      fstrDYNAMIC, fstrRESULT, fstrPARAM, &
      fstrCPL, fstrFREQ, hecLagMAT, restart_step_num)
    !C
    !C-- global variable
    !C
    type(hecmwst_local_mesh)             :: hecMESH
    type(hecmwst_matrix)                 :: hecMAT
    type(fstr_eigen)                     :: fstrEIG
    type(fstr_solid)                     :: fstrSOLID
    type(hecmwst_result_data)            :: fstrRESULT
    type(fstr_param)                     :: fstrPARAM
    type(fstr_dynamic)                   :: fstrDYNAMIC
    type(fstr_couple)                    :: fstrCPL
    type(fstr_freqanalysis)              :: fstrFREQ
    type(hecmwst_matrix_lagrange)        :: hecLagMAT
    integer(kind=kint)                   :: restart_step_num
 
    !C
    !C-- local variable
    !C
    integer(kind=kint)            :: numnode, numelm, startmode, endmode, nummode, ndof, im, in, ntotal, vistype
    integer(kind=kint)            :: numfreq, idnode, numdisp, imonit, nmonit, monitnode(1)
    integer(kind=kint)            :: freqiout(3)
    integer(kind=kint)            :: ierr
    integer(kind=kint), parameter :: ilogin = 9056
    real(kind=kreal), allocatable :: eigenvalue(:), loadvecre(:), loadvecim(:)
    real(kind=kreal), allocatable :: bjre(:), bjim(:)
    real(kind=kreal), allocatable :: dispre(:), dispim(:), velre(:), velim(:), accre(:), accim(:)
    real(kind=kreal)              :: freq, omega, val, dx, dy, dz, f_start, f_end
    real(kind=kreal)              :: t_start, t_end, time, dxi, dyi, dzi
    type(fstr_freqanalysis_data)  :: freqdata
 
    numnode   = hecmesh%n_node
    numelm    = hecmesh%n_elem
    ndof      = hecmesh%n_dof
    startmode = fstrfreq%start_mode
    endmode   = fstrfreq%end_mode
    nummode   = endmode - startmode +1
 
    freqdata%numMode    = nummode
    freqdata%numNodeDOF = numnode*ndof
    allocate(freqdata%eigOmega(nummode))
    allocate(freqdata%eigVector(numnode*ndof, nummode))
 
    call setupfreqparam(fstrdynamic, f_start, f_end, numfreq, freqdata%rayAlpha, freqdata%rayBeta, idnode, vistype, freqiout)
    call setupoutputitem(fstrsolid, freqiout)
 
    monitnode(1) = idnode
    imonit       = 0
    if( node_global_to_local(hecmesh, monitnode, 1) == 1 ) then
      if( monitnode(1) <= hecmesh%nn_internal ) imonit = monitnode(1)  !monitor is reported by the rank owning it
    end if
    nmonit = 0
    if( imonit > 0 ) nmonit = 1
    call hecmw_allreduce_i1(hecmesh, nmonit, hecmw_sum)
    if( nmonit == 0 ) then
      if( myrank == 0 ) then
        write(*,   *) "Error: monitor node not found:", idnode
        write(ilog,*) "Error: monitor node not found:", idnode
      end if
      call hecmw_abort( hecmw_comm_get_comm() )
    end if
 
    if( myrank == 0 ) then
      write(*,*) "Rayleigh alpha:", freqdata%rayAlpha
      write(*,*) "Rayleigh beta:", freqdata%rayBeta
      write(ilog,*) "Rayleigh alpha:", freqdata%rayAlpha
      write(ilog,*) "Rayleigh beta:", freqdata%rayBeta
    end if
 
 
    allocate(eigenvalue(nummode))
    allocate(loadvecre(numnode*ndof))
    allocate(loadvecim(numnode*ndof))
    allocate(bjre(nummode))
    allocate(bjim(nummode))
    allocate(dispre(numnode*ndof))
    allocate(dispim(numnode*ndof))
    allocate(velre(numnode*ndof))
    allocate(velim(numnode*ndof))
    allocate(accre(numnode*ndof))
    allocate(accim(numnode*ndof))
 
    loadvecre(:) = 0.0d0
    loadvecim(:) = 0.0d0
 
    if( myrank == 0 ) then
      write(*,*) "--frequency analysis--"
      write(*,   *) "read from=", trim(fstrfreq%eigenlog_filename)
      write(ilog,*) "read from=", trim(fstrfreq%eigenlog_filename)
      write(*,   *) "start mode=", startmode
      write(ilog,*) "start mode=", startmode
      write(*,   *) "end mode=", endmode
      write(ilog,*) "end mode=", endmode
    end if
    open(unit=ilogin, file=trim(fstrfreq%eigenlog_filename), status="OLD", action="READ", iostat=ierr)
    if( ierr /= 0 ) then
      if( myrank == 0 ) write(*,*) "Error: cannot open eigenlog file: ", trim(fstrfreq%eigenlog_filename)
      call hecmw_abort( hecmw_comm_get_comm() )
    endif
    call read_eigen_values(ilogin, startmode, endmode, eigenvalue, freqdata%eigOmega)
    !call read_eigen_vector(ilogin, startmode, endmode, ndof, numnode, freqData%eigVector)
    close(ilogin)
 
    call read_eigen_vector_res(hecmesh, startmode, endmode, ndof, numnode, freqdata%eigVector)
 
    call extract_surf2node(hecmesh, fstrfreq, ndof, loadvecre, loadvecim)
    call assemble_nodeload(hecmesh, fstrfreq, ndof, loadvecre, loadvecim)
 
    if( myrank == 0 ) write(*,*) "calc mass matrix"
    call calcmassmatrix(fstrparam, hecmesh, hecmat, fstrsolid, fstreig, heclagmat)
    if( myrank == 0 ) write(*,*) "scale eigenvector"
    call scaleeigenvector(hecmesh, fstreig, ndof*numnode, nummode, freqdata%eigVector)
 
    if( myrank == 0 ) then
      write(*,   *) "start frequency:", f_start
      write(ilog,*) "start frequency:", f_start
      write(*,   *) "end frequency:", f_end
      write(ilog,*) "end frequency:", f_end
      write(*   ,*) "number of the sampling points", numfreq
      write(ilog,*) "number of the sampling points", numfreq
      write(*   ,*) "monitor nodeid=", idnode
      write(ilog,*) "monitor nodeid=", idnode
    end if
 
    do im=1, numfreq
      freq = (f_end-f_start)/dble(numfreq)*dble(im) + f_start
      omega = 2.0d0 * 3.14159265358979d0 * freq
 
      call calcfreqcoeff(hecmesh, freqdata, loadvecre, loadvecim, omega, bjre, bjim)
      call calcdispvector(freqdata, bjre, bjim, dispre, dispim)
      call calcvelvector(freqdata, omega, bjre, bjim, velre, velim)
      call calcaccvector(freqdata, omega, bjre, bjim, accre, accim)
 
      if( imonit > 0 ) then
        dx  = sqrt(dispre(3*(imonit-1)+1)**2 + dispim(3*(imonit-1)+1)**2)
        dy  = sqrt(dispre(3*(imonit-1)+2)**2 + dispim(3*(imonit-1)+2)**2)
        dz  = sqrt(dispre(3*(imonit-1)+3)**2 + dispim(3*(imonit-1)+3)**2)
        val = sqrt(dx**2 + dy**2 + dz**2)
        write(*,    *) freq, "[Hz] : ", val
        write(ilog, *) freq, "[Hz] : ", val
      end if
 
      if( myrank == 0 ) then
        if(iresult==1) then
          write(*,   *) freq, "[Hz] : ", im, ".res"
          write(ilog,*) freq, "[Hz] : ", im, ".res"
        end if
        if(ivisual==1 .and. vistype==1) then
          write(*,   *) freq, "[Hz] : ", im, ".vis"
          write(ilog,*) freq, "[Hz] : ", im, ".vis"
        end if
      end if
      call output_result(hecmesh, hecmat, fstrsolid, fstrdynamic, 'fstrRES', 'frequency', freq, im, &
        iresult==1, ivisual==1 .and. vistype==1, dispre, dispim, velre, velim, accre, accim)
    end do
 
    call setupdynaparam(fstrdynamic, t_start, t_end, freq, numdisp)
    if( myrank == 0 ) then
      write(*,   *) "start time:", t_start
      write(ilog,*) "start time:", t_start
      write(*,   *) "end time:", t_end
      write(ilog,*) "end time:", t_end
      write(*,   *) "frequency:", freq
      write(ilog,*) "frequency:", freq
      write(*,   *) "node id:", idnode
      write(ilog,*) "node id:", idnode
      write(*,   *) "num disp:", numdisp
      write(ilog,*) "num disp:", numdisp
    end if
 
    omega = 2.0d0 * 3.14159265358979d0 * freq
    call calcfreqcoeff(hecmesh, freqdata, loadvecre, loadvecim, omega, bjre, bjim)
 
    do im=1, numdisp
      time = (t_end-t_start)/dble(numdisp)*dble(im-1) + t_start
      call calcdispvectortime(freqdata, time, omega, bjre, bjim, dispre, dispim)
      call calcvelvectortime(freqdata, time, omega, bjre, bjim, velre, velim)
      call calcaccvectortime(freqdata, time, omega, bjre, bjim, accre, accim)
      if( imonit > 0 ) then
        dx  = dispre(3*(imonit-1)+1)
        dy  = dispre(3*(imonit-1)+2)
        dz  = dispre(3*(imonit-1)+3)
        dxi = dispim(3*(imonit-1)+1)
        dyi = dispim(3*(imonit-1)+2)
        dzi = dispim(3*(imonit-1)+3)
      end if
 
      if( myrank == 0 ) then
        if(iresult==1) then
          write(*,   *) "time=", time, " : ", im, ".res"
          write(ilog,*) "time=", time, " : ", im, ".res"
        end if
        if(ivisual==1 .and. vistype==2) then
          write(*,   *) "time=", time, " : ", im, ".vis"
          write(ilog,*) "time=", time, " : ", im, ".vis"
        end if
      end if
      call output_result(hecmesh, hecmat, fstrsolid, fstrdynamic, 'fstrDYNA', 'TOTALTIME', time, im, &
        iresult==1, ivisual==1 .and. vistype==2, dispre, dispim, velre, velim, accre, accim)
    end do
 
    deallocate(freqdata%eigOmega)
    deallocate(freqdata%eigVector)
    deallocate(eigenvalue)
    deallocate(loadvecre)
    deallocate(loadvecim)
    deallocate(bjre)
    deallocate(bjim)
    deallocate(dispre)
    deallocate(dispim)
    deallocate(velre)
    deallocate(velim)
    deallocate(accre)
    deallocate(accim)
 
  end subroutine
 
  subroutine read_eigen_values(logfile, startmode, endmode, eigenvalue, anglfreq)
    !---- args
    integer(kind=kint), intent(in)  :: logfile
    integer(kind=kint), intent(in)  :: startmode
    integer(kind=kint), intent(in)  :: endmode
    real(kind=kreal), intent(inout) :: eigenvalue(:) !intend(endmode-startmode+1)
    real(kind=kreal), intent(inout) :: anglfreq(:)
    !---- vals
    integer(kind=kint)           :: im, endflag, id
    character(len=HECMW_MSG_LEN) :: line
    real(kind=kreal)             :: freq
    !---- body
 
    rewind(logfile)
    endflag = 0
    !Find eigenvalue header.
    do
      read(logfile, '(A80)', err=119) line
      if(trim(adjustl(line)) == "NO.  EIGENVALUE  FREQUENCY   (HZ)        X           Y           Z           X") then
        endflag = 1
        exit
      end if
    end do
    read(logfile, '(A80)') line
    !read eigenvalue
    do im=1, startmode-1
      read(logfile, '(A80)') line
    end do
    do im=1, (endmode-startmode+1)
      read(logfile, '(i5,3e12.4,a)', err=119) id, eigenvalue(im), anglfreq(im), freq, line
    end do
    return
 
    !error handling
    119 write(*,*) "Error to find eigenvalue information from logfile"
    write(ilog,*) "Error to find eigenvalue information from logfile"
    stop
  end subroutine
 
  subroutine read_eigen_vector(logfile, startmode, endmode, numdof, numnode, eigenvector)
    !---- args
    integer(kind=kint), intent(in)  :: logfile
    integer(kind=kint), intent(in)  :: startmode
    integer(kind=kint), intent(in)  :: endmode
    integer(kind=kint), intent(in)  :: numdof
    integer(kind=kint), intent(in)  :: numnode
    real(kind=kreal), intent(inout) :: eigenvector(:, :) !intend (numdof*NN,nmode)
    !---- vals
    integer(kind=kint)           :: im, in, gblid, j, idx
    real(kind=kreal)             :: vec(6)
    character(len=HECMW_MSG_LEN) :: line
    !---- body
 
    rewind(logfile)
    !Find first eigenvector header
    do im=1, startmode-1
      do
        read(logfile, '(a80)', err=119, end=119) line
        if(line(1:9) == " Mode No.") then
          exit ! goto next mode
        end if
      end do
    end do
 
    !read eigenvector
    do im=1, (endmode-startmode+1)
      !find header
      do
        read(logfile, '(a80)', err=119, end=119) line
        if(line(1:9) == " Mode No.") then
          exit !find eigenmode
        end if
      end do
 
      read(logfile, '(a80)', err=119) line
      read(logfile, '(a80)', err=119) line
      read(logfile, '(a80)', err=119) line
      read(logfile, '(a80)', err=119) line
      !read eigenvector component
      do in=1, numnode
        select case(numdof)
          case(2)
            read(logfile, '(i10,2e12.4)', err=119) gblid, (vec(j), j=1,2)
 
          case(3)
            !read(logfile, '(i10,3e12.4)', ERR=119) gblid, (vec(j), j=1,3)
            read(logfile, '(i10,3e16.8)', err=119) gblid, (vec(j), j=1,3)
 
          case(6)
            read(logfile, '(i10,6e12.4)', err=119) gblid, (vec(j), j=1,6)
          case default
            !error
            goto 119
        end select
 
        do j=1, numdof
          idx = (in-1)*numdof + j
          eigenvector(idx,im) = vec(j)
        end do
      end do
    end do
    return
 
    !error handling
    119 write(*,*) "Error to find eigenvector from logfile"
    write(ilog,*) "Error to find eigenvector from logfile"
    stop
 
  end subroutine
 
  subroutine read_eigen_vector_res(hecMESH, startmode, endmode, numdof, numnode, eigenvector)
    !---- args
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    integer(kind=kint), intent(in)       :: startmode
    integer(kind=kint), intent(in)       :: endmode
    integer(kind=kint), intent(in)       :: numdof
    integer(kind=kint), intent(in)       :: numnode
    real(kind=kreal), intent(inout)      :: eigenvector(:, :) !intend (numdof*NN,nmode)
    !---- vals
    integer(kind=kint), parameter :: compidx = 1 !Component index of displacement
    integer(kind=kint)            :: imode, idx, ind, a, b, nallcomp, j
    type(hecmwst_result_data)     :: eigenres
    character(len=HECMW_NAME_LEN) :: name
    !---- body
 
    name = 'result-in'
    do imode=startmode, endmode
      call nullify_result_data(eigenres)
      call hecmw_result_read_by_name(hecmesh, name, imode, eigenres)
 
      nallcomp = 0
      do ind=1,eigenres%nn_component
        nallcomp = nallcomp + eigenres%nn_dof(ind)
      end do
 
      idx = imode - startmode + 1
      do ind=1, numnode
        do j=1, numdof
          a = (ind-1)*nallcomp + j  !src vector index
          b = (ind-1)*numdof   + j
          eigenvector(b,idx) = eigenres%node_val_item(a)
        end do
      end do
      call free_result_data(eigenres)
    end do
 
  contains
 
    subroutine free_result_data(res)
      !---- args
      type(hecmwst_result_data) :: res
      !---- vals
      !---- body
      if(associated(res%nn_dof))        deallocate(res%nn_dof)
      if(associated(res%ne_dof))        deallocate(res%ne_dof)
      if(associated(res%node_label))    deallocate(res%node_label)
      if(associated(res%elem_label))    deallocate(res%elem_label)
      if(associated(res%node_val_item)) deallocate(res%node_val_item)
      if(associated(res%elem_val_item))  deallocate(res%elem_val_item)
    end subroutine
 
    subroutine nullify_result_data(res)
      !---- args
      type(hecmwst_result_data) :: res
      !---- vals
      !---- body
      nullify(res%nn_dof)
      nullify(res%ne_dof)
      nullify(res%node_label)
      nullify(res%elem_label)
      nullify(res%node_val_item)
      nullify(res%elem_val_item)
    end subroutine
 
 
  end subroutine
 
  !
  ! The displacement, the velocity and the acceleration of the frequency response
  ! analysis are selected by iout_list of !DYNAMIC, and the other items follow
  ! !OUTPUT_RES as in the static and the transient dynamic analysis
  !
  subroutine setupoutputitem(fstrSOLID, iout)
    !---- args
    type(fstr_solid), intent(inout) :: fstrSOLID
    integer(kind=kint), intent(in)  :: iout(3)
    !---- vals
    integer(kind=kint) :: i
    !---- body
 
    do i=3, 4
      fstrsolid%output_ctrl(i)%outinfo%on(1)  = (iout(1) == 1)
      fstrsolid%output_ctrl(i)%outinfo%on(15) = (iout(2) == 1)
      fstrsolid%output_ctrl(i)%outinfo%on(16) = (iout(3) == 1)
    end do
  end subroutine
 
  !
  ! Set one phase of the response to fstrSOLID and fstrDYNAMIC, and recover the
  ! strain and the stress of it with the routines of the static analysis.
  ! The incremental form keeps the incompatible mode of the 361 element consistent
  ! with the displacement, as the material is linear elastic here.
  !
  subroutine updatesolidresult(hecMESH, hecMAT, fstrSOLID, fstrDYNAMIC, istp, disp, vel, acc)
    !---- args
    type(hecmwst_local_mesh)       :: hecMESH
    type(hecmwst_matrix)           :: hecMAT
    type(fstr_solid)               :: fstrSOLID
    type(fstr_dynamic)             :: fstrDYNAMIC
    integer(kind=kint), intent(in) :: istp
    real(kind=kreal), intent(in)   :: disp(:) !intend (numnodeDOF)
    real(kind=kreal), intent(in)   :: vel(:) !intend (numnodeDOF)
    real(kind=kreal), intent(in)   :: acc(:) !intend (numnodeDOF)
    !---- vals
    integer(kind=kint) :: ndof, ntotal, idx
    !---- body
 
    ndof   = hecmesh%n_dof
    ntotal = ndof*hecmesh%n_node
 
    idx = 1
    if( fstrdynamic%idx_eqa == 1 .and. istp > 0 ) idx = 2
    fstrdynamic%DISP(1:ntotal,idx) = disp(1:ntotal)
    fstrdynamic%VEL (1:ntotal,idx) = vel(1:ntotal)
    fstrdynamic%ACC (1:ntotal,idx) = acc(1:ntotal)
 
    hecmat%X(1:ntotal) = disp(1:ntotal) - fstrsolid%unode(1:ntotal) - fstrsolid%dunode(1:ntotal)
    call fstr_apply_solution_increment( hecmesh, fstrsolid, ndof, hecmat%X )
    call fstr_updatenewton( hecmesh, hecmat, fstrsolid, 0.0d0, 0.0d0, 1 )
    call fstr_update_reaction_spc( 1, hecmesh, fstrsolid )
 
    if( ndof == 2 ) then
      call fstr_nodalstress2d( hecmesh, fstrsolid )
    else if( ndof == 6 ) then
      call fstr_nodalstress6d( hecmesh, fstrsolid )
    else
      call fstr_nodalstress3d( hecmesh, fstrsolid )
    end if
  end subroutine
 
  !
  ! The response is complex, so the real part and the imaginary part are written as
  ! two sets of the same items, distinguished by the label suffix. Both are recovered
  ! once and shared by the result file and the visualization data.
  !
  subroutine output_result(hecMESH, hecMAT, fstrSOLID, fstrDYNAMIC, nameid, glabel, gvalue, istp, &
      wres, wvis, dispRe, dispIm, velRe, velIm, accRe, accIm)
    !---- args
    type(hecmwst_local_mesh)       :: hecMESH
    type(hecmwst_matrix)           :: hecMAT
    type(fstr_solid)               :: fstrSOLID
    type(fstr_dynamic)             :: fstrDYNAMIC
    character(len=*), intent(in)   :: nameid
    character(len=*), intent(in)   :: glabel
    real(kind=kreal), intent(in)   :: gvalue
    integer(kind=kint), intent(in) :: istp
    logical, intent(in)            :: wres, wvis
    real(kind=kreal), intent(in)   :: dispre(:), dispim(:) !intend (numnodeDOF)
    real(kind=kreal), intent(in)   :: velre(:), velim(:) !intend (numnodeDOF)
    real(kind=kreal), intent(in)   :: accre(:), accim(:) !intend (numnodeDOF)
    !---- vals
    type(hecmwst_result_data)       :: visre, visim
    character(len=HECMW_HEADER_LEN) :: header
    character(len=HECMW_MSG_LEN)    :: comment
    character(len=HECMW_NAME_LEN)   :: label, resid
    real(kind=kreal)                :: gval(1)
    !---- body
 
    if( .not. (wres .or. wvis) ) return
 
    if( wres ) then
      header='*fstrresult'
      comment='frequency_result'
      call hecmw_result_init(hecmesh, istp, header, comment)
 
      label   = glabel
      gval(1) = gvalue
      call hecmw_result_add(hecmw_result_dtype_global, 1, label, gval)
    end if
 
    call updatesolidresult(hecmesh, hecmat, fstrsolid, fstrdynamic, istp, dispre, velre, accre)
    if( wres ) call fstr_write_result_add(hecmesh, fstrsolid, istp, '_real', fstrdynamic)
    if( wvis ) call fstr_make_result(hecmesh, fstrsolid, visre, istp, gvalue, fstrdynamic, '_real')
 
    call updatesolidresult(hecmesh, hecmat, fstrsolid, fstrdynamic, istp, dispim, velim, accim)
    if( wres ) call fstr_write_result_add(hecmesh, fstrsolid, istp, '_imag', fstrdynamic)
    if( wvis ) call fstr_make_result(hecmesh, fstrsolid, visim, istp, gvalue, fstrdynamic, '_imag')
 
    if( wres ) then
      resid = nameid
      call hecmw_result_write_by_name(resid)
      call hecmw_result_finalize()
    end if
 
    if( wvis ) then
      call fstr_merge_result(hecmesh, visre, visim)
      call hecmw_result_free(visim)
      call fstr2hecmw_mesh_conv(hecmesh)
      call hecmw_visualize_init
      call hecmw_visualize( hecmesh, visre, istp )
      call hecmw2fstr_mesh_conv(hecmesh)
      call hecmw_result_free(visre)
    end if
    return
  end subroutine
 
  subroutine extract_surf2node(hecMESH, freqData, numdof, loadvecRe, loadvecIm)
    !---- args
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    type(fstr_freqanalysis), intent(in)  :: freqData
    integer(kind=kint), intent(in)       :: numdof
    real(kind=kreal), intent(inout)      :: loadvecre(:) !intend(numnode*ndof)
    real(kind=kreal), intent(inout)      :: loadvecim(:) !intend(numnode*ndof)
    !---- vals
    integer(kind=kint), parameter   :: MAXNODE = 100
    integer(kind=kint)              :: sgrpID, is, ie, ic, nsurf, ic_type, outtype, node_index(MAXNODE)
    integer(kind=kint)              :: nn, iss, nodeid, dof_index, ndof
    integer(kind=kint)              :: i, j, k, l, m, isn, nsize
    integer(kind=kint)              :: iwk(60), nodLOCAL(20)
    real(kind=kreal)                :: vect(60), xx(20), yy(20), zz(20), forcere(3), forceim(3)
    !---- body
 
    ndof = 3
    do i=1,freqdata%FLOAD_ngrp_tot
      if(freqdata%FLOAD_ngrp_TYPE(i) == kfloadtype_surf) then  !FLOAD type=surface
        sgrpid     = freqdata%FLOAD_ngrp_ID(i)
        dof_index  = freqdata%FLOAD_ngrp_DOF(i)
        forcere(:) = 0.0d0
        forceim(:) = 0.0d0
        forcere(dof_index) = freqdata%FLOAD_ngrp_valre(i)
        forceim(dof_index) = freqdata%FLOAD_ngrp_valim(i)
 
        is = hecmesh%surf_group%grp_index(sgrpid-1) + 1
        ie = hecmesh%surf_group%grp_index(sgrpid)
        do j=is, ie
          ic      = hecmesh%surf_group%grp_item(2*j-1)
          nsurf   = hecmesh%surf_group%grp_item(2*j)
          ic_type = hecmesh%elem_type(ic)
          nn  = hecmw_get_max_node(ic_type)
          isn = hecmesh%elem_node_index(ic-1)
          do k=1, nn
            nodlocal(k) = hecmesh%elem_node_item(isn+k)
            xx(k) = hecmesh%node(3*nodlocal(k)-2)
            yy(k) = hecmesh%node(3*nodlocal(k)-1)
            zz(k) = hecmesh%node(3*nodlocal(k)  )
            do l=1, ndof
              iwk(ndof*(k-1)+l) = ndof*(nodlocal(k)-1)+l
            end do
          end do
 
          call dl_c3_freq(ic_type, nn, xx, yy, zz, nsurf, forcere, vect, nsize)
          do k=1,nsize
            loadvecre(iwk(k)) = loadvecre(iwk(k)) + vect(k)
          end do
 
          call dl_c3_freq(ic_type, nn, xx, yy, zz, nsurf, forceim, vect, nsize)
          do k=1,nsize
            loadvecim(iwk(k)) = loadvecim(iwk(k)) + vect(k)
          end do
        end do
      end if
    end do
 
    return
  end subroutine
 
  subroutine dl_c3_freq(ETYPE, NN, XX, YY, ZZ, LTYPE, force, VECT, nsize)
    !---- args
    integer(kind=kint), intent(in)    :: ETYPE !--solid element type
    integer(kind=kint), intent(in)    :: NN !--node num
    integer(kind=kint), intent(in)    :: LTYPE !--solid element face
    real(kind=kreal), intent(in)      :: xx(:) !--node x pos
    real(kind=kreal), intent(in)      :: yy(:) !--node y pos
    real(kind=kreal), intent(in)      :: zz(:) !--node z pos
    real(kind=kreal), intent(in)      :: force(3) !--node surfforce
    real(kind=kreal), intent(inout)   :: vect(:)
    integer(kind=kint), intent(inout) :: nsize
    !---- vals
    integer(kind=kint), parameter :: NDOF = 3
    real(kind=kreal)              :: wg, area
    integer(kind=kint)            :: NOD(NN)
    real(kind=kreal)              :: elecoord(3, nn), localcoord(3), normal(3)
    real(kind=kreal)              :: h(nn)
    integer(kind=kint)            :: I, IG2, NSUR, SURTYPE
    !---- body
 
    call getsubface( etype, ltype, surtype, nod )
    nsur = getnumberofnodes( surtype )
 
    do i=1,nsur
      elecoord(1,i)=xx(nod(i))
      elecoord(2,i)=yy(nod(i))
      elecoord(3,i)=zz(nod(i))
    end do
    nsize         = nn*ndof
    vect(1:nsize) = 0.0d0
    do ig2=1,numofquadpoints( surtype )
      call getquadpoint( surtype, ig2, localcoord(1:2) )
      call getshapefunc( surtype, localcoord(1:2), h(1:nsur) )
 
      wg=getweight( surtype, ig2 )
      normal=surfacenormal( surtype, nsur, localcoord(1:2), elecoord(:,1:nsur) )
      area=dsqrt( dot_product(normal, normal) )
      do i=1,nsur
        vect(3*nod(i)-2)=vect(3*nod(i)-2)+area*wg*h(i)*force(1)
        vect(3*nod(i)-1)=vect(3*nod(i)-1)+area*wg*h(i)*force(2)
        vect(3*nod(i)  )=vect(3*nod(i)  )+area*wg*h(i)*force(3)
      end do
    end do
  end subroutine
 
  subroutine assemble_nodeload(hecMESH, freqData, numdof, loadvecRe, loadvecIm)
    !---- args
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    type(fstr_freqanalysis), intent(in)  :: freqData
    integer(kind=kint), intent(in)       :: numdof
    real(kind=kreal), intent(inout)      :: loadvecre(:)
    real(kind=kreal), intent(inout)      :: loadvecim(:)
    !---- vals
    integer(kind=kint) :: i, vecsize, ig, is, ie, in, nodeid, dof_index
 
    !---- body
 
    do i=1, freqdata%FLOAD_ngrp_tot
      if(freqdata%FLOAD_ngrp_TYPE(i) == kfloadtype_node) then
        ig = freqdata%FLOAD_ngrp_ID(i)
        is = hecmesh%node_group%grp_index(ig-1) + 1
        ie = hecmesh%node_group%grp_index(ig)
        do in=is, ie
          nodeid    = hecmesh%node_group%grp_item(in)
          dof_index = freqdata%FLOAD_ngrp_DOF(i)
          loadvecre((nodeid-1)*numdof + dof_index) = loadvecre((nodeid-1)*numdof + dof_index) + freqdata%FLOAD_ngrp_valre(i)
          loadvecim((nodeid-1)*numdof + dof_index) = loadvecim((nodeid-1)*numdof + dof_index) + freqdata%FLOAD_ngrp_valim(i)
        end do
      end if
    end do
    return
  end subroutine
 
  subroutine calcmassmatrix(fstrPARAM, hecMESH, hecMAT, fstrSOLID, fstrEIG, hecLagMAT)
    !---- args
    type(fstr_param), intent(in)         :: fstrPARAM
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    type(hecmwst_matrix), intent(inout)  :: hecMAT
    type(fstr_solid), intent(inout)      :: fstrSOLID
    type(fstr_eigen), intent(inout)      :: fstrEIG
    type(hecmwst_matrix_lagrange), intent(inout) :: hecLagMAT
    !---- vals
    integer(kind=kint) :: ntotal
    type(hecmwst_ebc) :: hecEBC
    !---- body
 
 
    fstrsolid%dunode = 0.d0
    call fstr_creatematrix_and_dampingforce( hecmesh, hecmat, fstrsolid, 0.d0, 0.d0 )
    call hecmw_ebc_init(hecmat, hecebc)
    call fstr_addbc(1, hecmesh, hecmat, fstrsolid, fstrparam, heclagmat, 2, hecebc)
    call hecmw_ebc_apply(hecmesh, hecmat, hecebc)
    call hecmw_ebc_finalize(hecebc)
 
    call setmass(fstrsolid, hecmesh, hecmat, fstreig)
 
  end subroutine
 
  subroutine scaleeigenvector(hecMESH, fstrEIG, ntotaldof, nmode, eigenvector)
    !---- args
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    type(fstr_eigen), intent(in)    :: fstrEIG
    integer(kind=kint), intent(in)  :: ntotaldof
    integer(kind=kint), intent(in)  :: nmode
    real(kind=kreal), intent(inout) :: eigenvector(:, :)
    !---- vals
    integer(kind=kint) :: imode, idof, nintdof
    real(kind=kreal)   :: mas
    !---- body
 
    nintdof = hecmesh%nn_internal*hecmesh%n_dof
 
    do imode=1,nmode
      mas = 0.0d0
      do idof=1,nintdof
        mas = mas + fstreig%mass(idof)*eigenvector(idof,imode)**2
      end do
      call hecmw_allreduce_r1(hecmesh, mas, hecmw_sum)
      do idof=1,ntotaldof
        eigenvector(idof,imode) =  eigenvector(idof,imode) / sqrt(mas)
      end do
    end do
  end subroutine
 
  subroutine checkorthvector(fstrEIG, eigenvector, imode, jmode, prod)
    !---- args
    type(fstr_eigen), intent(in)    :: fstreig
    real(kind=kreal), intent(in)    :: eigenvector(:, :)
    integer(kind=kint), intent(in)  :: imode
    integer(kind=kint), intent(in)  :: jmode
    real(kind=kreal), intent(inout) :: prod
    !---- vals
    integer(kind=kint) :: idof, s
    !---- body
    s = size(eigenvector(:,1))
    prod = 0.0d0
 
    do idof=1, s
      prod = prod + eigenvector(idof,imode)*fstreig%mass(idof)*eigenvector(idof,jmode)
    end do
    return
  end subroutine
 
  subroutine writeoutvector(im, vector)
    !---- args
    integer(kind=kint), intent(in) :: im
    real(kind=kreal), intent(in)   :: vector(:)
    !---- vals
    integer(kind=kint) :: i, s
    !---- body
    s = size(vector)
    do i=1, s
      if(i == 1) then
        write(*,'("eigenvec",i2.2,":[",e12.5,", ")') im, vector(i)
      else if(i /= s) then
        write(*,'(e12.5,", ")') vector(i)
      else
        write(*,'(e12.5,"];")') vector(i)
      end if
    end do
    write(*,*)
    return
  end subroutine
 
  subroutine calcdotproduct(hecMESH, a, b, c)
    !---- args
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    real(kind=kreal), intent(in)    :: a(:)
    real(kind=kreal), intent(in)    :: b(:)
    real(kind=kreal), intent(inout) :: c
    !---- vals
    integer(kind=kint) :: nintdof
    !---- body
    nintdof = hecmesh%nn_internal*hecmesh%n_dof  !external nodes are owned by another rank
    c = dot_product(a(1:nintdof), b(1:nintdof))
    call hecmw_allreduce_r1(hecmesh, c, hecmw_sum)
    return
  end subroutine
 
  subroutine calcfreqcoeff(hecMESH, freqData, loadRe, loadIm, inpOmega, bjRe, bjIm)
    !---- args
    type(hecmwst_local_mesh), intent(in)     :: hecMESH
    type(fstr_freqanalysis_data), intent(in) :: freqData
    real(kind=kreal), intent(in)             :: loadre(:) !intend (numNodeDOF)
    real(kind=kreal), intent(in)             :: loadim(:) !intend (numNodeDOF)
    real(kind=kreal), intent(in)             :: inpomega
    real(kind=kreal), intent(inout)          :: bjre(:) !intend (numMode)
    real(kind=kreal), intent(inout)          :: bjim(:) !intend (numMode)
    !---- vals
    integer(kind=kint) :: imode
    real(kind=kreal)   :: ujfr, ujfi, a, b, alp, beta
    !---- body
 
    alp  = freqdata%rayAlpha
    beta = freqdata%rayBeta
 
    do imode=1, freqdata%numMode
      call calcdotproduct(hecmesh, freqdata%eigVector(:,imode), loadre, ujfr)
      call calcdotproduct(hecmesh, freqdata%eigVector(:,imode), loadim, ujfi)
 
      a = ujfr*(freqdata%eigOmega(imode)**2 - inpomega**2) + ujfi*(alp + beta*freqdata%eigOmega(imode)**2)*inpomega
      b = (freqdata%eigOmega(imode)**2 - inpomega**2)**2 + ((alp + beta*freqdata%eigOmega(imode)**2)*inpomega)**2
      bjre(imode) = a / b
 
      a = ujfi*(freqdata%eigOmega(imode)**2 -inpomega**2) - ujfr*(alp + beta*freqdata%eigOmega(imode)**2)*inpomega
      b = (freqdata%eigOmega(imode)**2 - inpomega**2)**2 + ((alp + beta*freqdata%eigOmega(imode)**2)*inpomega)**2
      bjim(imode) = a / b
    end do
    return
  end subroutine
 
  subroutine calcdispvector(freqData, bjRe, bjIm, dispRe, dispIm)
    !---- args
    type(fstr_freqanalysis_data), intent(in) :: freqdata
    real(kind=kreal), intent(in)             :: bjre(:) !intend (numMode)
    real(kind=kreal), intent(in)             :: bjim(:) !intend (numMode)
    real(kind=kreal), intent(inout)          :: dispre(:) !intend (numNodeDOF)
    real(kind=kreal), intent(inout)          :: dispim(:) !intend (numNodeDOF)
    !---- vals
    integer(kind=kint) :: imode
    !---- body
 
    dispre(:) = 0.0d0
    dispim(:) = 0.0d0
 
    do imode=1, freqdata%numMode
      dispre(:) = dispre(:) + bjre(imode)*freqdata%eigVector(:,imode)
      dispim(:) = dispim(:) + bjim(imode)*freqdata%eigVector(:,imode)
    end do
    return
  end subroutine
 
  subroutine calcvelvector(freqData, omega, bjRe, bjIm, velRe, velIm)
    !---- args
    type(fstr_freqanalysis_data), intent(in) :: freqData
    real(kind=kreal), intent(in)             :: omega
    real(kind=kreal), intent(in)             :: bjre(:) !intend (numMode)
    real(kind=kreal), intent(in)             :: bjim(:) !intend (numMode)
    real(kind=kreal), intent(inout)          :: velre(:) !intend (numNodeDOF)
    real(kind=kreal), intent(inout)          :: velim(:) !intend (numNodeDOF)
    !---- vals
    integer(kind=kint) :: imode
    !---- body
 
    velre(:) = 0.0d0
    velim(:) = 0.0d0
 
    do imode=1, freqdata%numMode
      velre(:) = velre(:) - omega * bjim(imode) * freqdata%eigVector(:,imode)
      velim(:) = velim(:) + omega * bjre(imode) * freqdata%eigVector(:,imode)
    end do
  end subroutine
 
  subroutine calcaccvector(freqData, omega, bjRe, bjIm, accRe, accIm)
    !---- args
    type(fstr_freqanalysis_data), intent(in) :: freqData
    real(kind=kreal), intent(in)             :: omega
    real(kind=kreal), intent(in)             :: bjre(:) !intend (numMode)
    real(kind=kreal), intent(in)             :: bjim(:) !intend (numMode)
    real(kind=kreal), intent(inout)          :: accre(:) !intend (numNodeDOF)
    real(kind=kreal), intent(inout)          :: accim(:) !intend (numNodeDOF)
    !---- vals
    integer(kind=kint) :: imode
    !---- body
 
    accre(:) = 0.0d0
    accim(:) = 0.0d0
 
    do imode=1, freqdata%numMode
      accre(:) = accre(:) - omega**2 * bjre(imode) * freqdata%eigVector(:,imode)
      accim(:) = accim(:) - omega**2 * bjim(imode) * freqdata%eigVector(:,imode)
    end do
 
  end subroutine
 
  subroutine setupfreqparam(fstrDYNAMIC, f_start, f_end, numfreq, raym, rayk, idnode, vistype, ioutl)
    !---- args
    type(fstr_dynamic), intent(in)    :: fstrDYNAMIC
    real(kind=kreal), intent(inout)   :: f_start
    real(kind=kreal), intent(inout)   :: f_end
    integer(kind=kint), intent(inout) :: numfreq
    real(kind=kreal), intent(inout)   :: raym
    real(kind=kreal), intent(inout)   :: rayk
    integer(kind=kint), intent(inout) :: idnode
    integer(kind=kint), intent(inout) :: vistype
    integer(kind=kint), intent(inout) :: ioutl(3)
    !---- vals
 
    !---- body
    f_start    = fstrdynamic%t_start
    f_end      = fstrdynamic%t_end
    numfreq    = fstrdynamic%n_step
    raym       = fstrdynamic%ray_m
    rayk       = fstrdynamic%ray_k
    idnode     = fstrdynamic%nout_monit
    vistype    = fstrdynamic%ngrp_monit
    ioutl(1:3) = fstrdynamic%iout_list(1:3)
    return
  end subroutine
 
  subroutine calcdispvectortime(freqData, time, omega, bjRe, bjIm, dispRe, dispIm)
    !---- args
    type(fstr_freqanalysis_data), intent(in) :: freqData
    real(kind=kreal), intent(in)             :: time
    real(kind=kreal), intent(in)             :: omega
    real(kind=kreal), intent(in)             :: bjre(:) !intend (numMode)
    real(kind=kreal), intent(in)             :: bjim(:) !intend (numMode)
    real(kind=kreal), intent(inout)          :: dispre(:) !intend (numNodeDOF)
    real(kind=kreal), intent(inout)          :: dispim(:) !intend (numNodeDOF)
    !---- vals
    integer(kind=kint)  :: imode, idf, s
    complex(kind=kreal) :: a, b, c
    !---- body
 
    dispre(:) = 0.0d0
    dispim(:) = 0.0d0
    a = exp(dcmplx(0.0d0, omega*time))
 
    do imode=1, freqdata%numMode
      s = size(freqdata%eigvector(:,imode))
      b = dcmplx(bjre(imode), bjim(imode)) * a
      do idf=1, s
        c = b*dcmplx(freqdata%eigVector(idf,imode), 0.0d0)
        dispre(idf) = dispre(idf) + dble(c)
        dispim(idf) = dispim(idf) + imag(c)
      end do
    end do
    return
  end subroutine
 
  subroutine calcvelvectortime(freqData, time, omega, bjRe, bjIm, velRe, velIm)
    !---- args
    type(fstr_freqanalysis_data), intent(in) :: freqData
    real(kind=kreal), intent(in)             :: time
    real(kind=kreal), intent(in)             :: omega
    real(kind=kreal), intent(in)             :: bjre(:) !intend (numMode)
    real(kind=kreal), intent(in)             :: bjim(:) !intend (numMode)
    real(kind=kreal), intent(inout)          :: velre(:) !intend (numNodeDOF)
    real(kind=kreal), intent(inout)          :: velim(:) !intend (numNodeDOF)
    !---- vals
    integer(kind=kint)  :: imode, idf, s
    complex(kind=kreal) :: a, b, c
    !---- body
 
    velre(:) = 0.0d0
    velim(:) = 0.0d0
    a = dcmplx(0.0d0, 1.0d0)*dcmplx(omega, 0.0d0)*exp(dcmplx(0.0d0, omega*time))
 
    do imode=1, freqdata%numMode
      s = size(freqdata%eigvector(:,imode))
      b = dcmplx(bjre(imode), bjim(imode)) * a
      do idf=1, s
        c = b*dcmplx(freqdata%eigVector(idf,imode), 0.0d0)
        velre(idf) = velre(idf) + dble(c)
        velim(idf) = velim(idf) + imag(c)
      end do
    end do
    return
  end subroutine
 
  subroutine calcaccvectortime(freqData, time, omega, bjRe, bjIm, accRe, accIm)
    !---- args
    type(fstr_freqanalysis_data), intent(in) :: freqData
    real(kind=kreal), intent(in)             :: time
    real(kind=kreal), intent(in)             :: omega
    real(kind=kreal), intent(in)             :: bjre(:) !intend (numMode)
    real(kind=kreal), intent(in)             :: bjim(:) !intend (numMode)
    real(kind=kreal), intent(inout)          :: accre(:) !intend (numNodeDOF)
    real(kind=kreal), intent(inout)          :: accim(:) !intend (numNodeDOF)
    !---- vals
    integer(kind=kint)  :: imode, idf, s
    complex(kind=kreal) :: a, b, c
    !---- body
 
    accre(:) = 0.0d0
    accim(:) = 0.0d0
    a = dcmplx(-1.0d0, 0.0d0)*dcmplx(omega**2, 0.0d0)*exp(dcmplx(0.0d0, omega*time))
 
    do imode=1, freqdata%numMode
      s = size(freqdata%eigvector(:,imode))
      b = dcmplx(bjre(imode), bjim(imode)) * a
      do idf=1, s
        c = b*dcmplx(freqdata%eigVector(idf,imode), 0.0d0)
        accre(idf) = accre(idf) + dble(c)
        accim(idf) = accim(idf) + imag(c)
      end do
    end do
    return
  end subroutine
 
  subroutine setupdynaparam(fstrDYNAMIC, t_start, t_end, dynafreq, numdisp)
    !---- args
    type(fstr_dynamic), intent(in)    :: fstrDYNAMIC
    real(kind=kreal), intent(inout)   :: t_start
    real(kind=kreal), intent(inout)   :: t_end
    real(kind=kreal), intent(inout)   :: dynafreq
    integer(kind=kint), intent(inout) :: numdisp
    !---- vals
    !---- body
    t_start  = fstrdynamic%gamma
    t_end    = fstrdynamic%beta
    dynafreq = fstrdynamic%t_delta
    numdisp  = fstrdynamic%nout
    return
  end subroutine
 
end module