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_visout
  use m_fstr
 
  implicit none
contains
  subroutine fstr_freq_result_init(hecMESH, numcomp, fstrRESULT)
    !---- args
    type(hecmwst_local_mesh), intent(in)     :: hecMESH
    integer(kind=kint), intent(in)           :: numcomp
    type(hecmwst_result_data), intent(inout) :: fstrRESULT
    !---- vals
    !---- body
 
    call hecmw_nullify_result_data(fstrresult)
    fstrresult%ng_component = 0
    fstrresult%nn_component = numcomp
    fstrresult%ne_component = 0
    allocate( fstrresult%nn_dof(numcomp) )
    allocate( fstrresult%node_label(numcomp) )
    allocate( fstrresult%node_val_item(numcomp*hecmesh%n_dof*hecmesh%n_node))  ! Should we use nn_internal?
  end subroutine
 
  subroutine fstr_freq_result_add(fstrRESULT, hecMESH, comp_index, ndof, label, vect)
    !---- args
    type(hecmwst_result_data), intent(inout)  :: fstrRESULT
    type(hecmwst_local_mesh), intent(in)      :: hecMESH
    integer(kind=kint), intent(in)            :: comp_index
    integer(kind=kint), intent(in)            :: ndof
    character(len=HECMW_NAME_LEN), intent(in) :: label
    real(kind=kreal), intent(in)              :: vect(:)
    !---- vals
    integer(kind=kint) :: i, k, alldof, offset
    !---- body
 
    fstrresult%nn_dof(comp_index)     = ndof
    fstrresult%node_label(comp_index) = label
    alldof = fstrresult%nn_component*ndof
    offset = ndof*(comp_index-1)
    do i=1, hecmesh%n_node
      do k=1, ndof
        fstrresult%node_val_item(alldof*(i-1) + k + offset) = vect(ndof*(i-1) + k)
      end do
    end do
 
  end subroutine
 
end module
 
 
module fstr_frequency_analysis
 
  use m_fstr
  use m_fstr_stiffmatrix
  use m_fstr_addbc
  use m_fstr_eig_setmass
  use fstr_frequency_visout
  use m_hecmw2fstr_mesh_conv
 
  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
    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(:), dvare(:), dvaim(:), disp(:), vel(:), acc(:)
    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%nn_internal
    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)
    write(*,*) "Rayleigh alpha:", freqdata%rayAlpha
    write(*,*) "Rayleigh beta:", freqdata%rayBeta
    write(ilog,*) "Rayleigh alpha:", freqdata%rayAlpha
    write(ilog,*) "Rayleigh beta:", freqdata%rayBeta
 
 
    allocate(eigenvalue(nummode))
    allocate(loadvecre(numnode*ndof))
    allocate(loadvecim(numnode*ndof))
    allocate(bjre(nummode))
    allocate(bjim(nummode))
    allocate(dvare(numnode*ndof))
    allocate(dvaim(numnode*ndof))
    allocate(disp(numnode*ndof))
    allocate(vel(numnode*ndof))
    allocate(acc(numnode*ndof))
    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
 
    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
    open(unit=ilogin, file=trim(fstrfreq%eigenlog_filename), status="OLD", action="READ", iostat=ierr)
    if( ierr /= 0 ) then
      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)
 
    write(*,*) "calc mass matrix"
    call calcmassmatrix(fstrparam, hecmesh, hecmat, fstrsolid, fstreig, heclagmat)
    write(*,*) "scale eigenvector"
    call scaleeigenvector(fstreig, ndof*numnode, nummode, freqdata%eigVector)
 
    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
 
    do im=1, numfreq
      freq = (f_end-f_start)/dble(numfreq)*dble(im) + f_start
      omega = 2.0d0 * 3.14159265358979d0 * freq
 
      call calcfreqcoeff(freqdata, loadvecre, loadvecim, omega, bjre, bjim)
      call calcdispvector(freqdata, bjre, bjim, dvare, dvaim)
 
      dx  = sqrt(dvare(3*(idnode-1)+1)**2 + dvaim(3*(idnode-1)+1)**2)
      dy  = sqrt(dvare(3*(idnode-1)+2)**2 + dvaim(3*(idnode-1)+2)**2)
      dz  = sqrt(dvare(3*(idnode-1)+3)**2 + dvaim(3*(idnode-1)+3)**2)
      val = sqrt(dx**2 + dy**2 + dz**2)
      write(*,    *) freq, "[Hz] : ", val
      write(ilog, *) freq, "[Hz] : ", val
      disp(:) = abs(cmplx(dvare(:), dvaim(:)))
 
      call calcvelvector(freqdata, omega, bjre, bjim, dvare, dvaim)
      vel(:) = abs(cmplx(dvare(:), dvaim(:)))
 
      call calcaccvector(freqdata, omega, bjre, bjim, dvare, dvaim)
      acc(:) = abs(cmplx(dvare(:), dvaim(:)))
 
      if(iresult==1) then
        write(*,   *) freq, "[Hz] : ", im, ".res"
        write(ilog,*) freq, "[Hz] : ", im, ".res"
        call output_resfile(hecmesh, freq, im, disp, vel, acc, freqiout)
      end if
      if(ivisual==1 .and. vistype==1) then
        write(*,   *) freq, "[Hz] : ", im, ".vis"
        write(ilog,*) freq, "[Hz] : ", im, ".vis"
        call output_visfile(hecmesh, im, disp, vel, acc, freqiout)
      end if
    end do
 
    call setupdynaparam(fstrdynamic, t_start, t_end, freq, numdisp)
    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
 
    omega = 2.0d0 * 3.14159265358979d0 * freq
    call calcfreqcoeff(freqdata, loadvecre, loadvecim, omega, bjre, bjim)
    call calcdispvector(freqdata, bjre, bjim, dvare, dvaim)
 
    do im=1, numdisp
      time = (t_end-t_start)/dble(numdisp)*dble(im-1) + t_start
      call calcdispvectortime(freqdata, time, omega, bjre, bjim, dvare, dvaim)
      dx  = dvare(3*(idnode-1)+1)
      dy  = dvare(3*(idnode-1)+2)
      dz  = dvare(3*(idnode-1)+3)
      dxi = dvaim(3*(idnode-1)+1)
      dyi = dvaim(3*(idnode-1)+2)
      dzi = dvaim(3*(idnode-1)+3)
 
      call calcvelvectortime(freqdata, time, omega, bjre, bjim, velre, velim)
      call calcaccvectortime(freqdata, time, omega, bjre, bjim, accre, accim)
      if(iresult==1) then
        write(*,   *) "time=", time, " : ", im, ".res"
        write(ilog,*) "time=", time, " : ", im, ".res"
        call outputdyna_resfile(hecmesh, time, im, dvare, dvaim, velre, velim, accre, accim, freqiout)
      end if
      if(ivisual==1 .and. vistype==2) then
        write(*,   *) "time=", time, " : ", im, ".vis"
        write(ilog,*) "time=", time, " : ", im, ".vis"
        call outputdyna_visfile(hecmesh, im, dvare, dvaim, velre, velim, accre, accim, freqiout)
      end if
    end do
 
    deallocate(freqdata%eigOmega)
    deallocate(freqdata%eigVector)
    deallocate(eigenvalue)
    deallocate(loadvecre)
    deallocate(loadvecim)
    deallocate(bjre)
    deallocate(bjim)
    deallocate(dvare)
    deallocate(dvaim)
    deallocate(disp)
    deallocate(vel)
    deallocate(acc)
    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,imode) = 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
 
  subroutine output_resfile(hecMESH, freq, ifreq, disp, vel, acc, iout)
    !---- args
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    real(kind=kreal), intent(in)         :: freq
    integer(kind=kint), intent(in)       :: ifreq
    real(kind=kreal), intent(in)         :: disp(:) !intend (numnodeDOF)
    real(kind=kreal), intent(in)         :: vel(:) !intend (numnodeDOF)
    real(kind=kreal), intent(in)         :: acc(:) !intend (numnodeDOF)
    integer(kind=kint), intent(in)       :: iout(3)
    !---- vals
    integer(kind=kint)              :: im
    character(len=HECMW_HEADER_LEN) :: header
    character(len=HECMW_MSG_LEN)    :: comment
    character(len=HECMW_NAME_LEN)   :: label, nameid
    real(kind=kreal)                :: freqval(1)
    !---- body
 
    nameid='fstrRES'
    header='*fstrresult'
    comment='frequency_result'
    call hecmw_result_init(hecmesh, ifreq, header, comment)
 
    label = "frequency"
    freqval(1) = freq
    call hecmw_result_add(hecmw_result_dtype_global, 1, label, freqval)
 
    if(iout(1) == 1) then
      label='displacement'
      call hecmw_result_add(hecmw_result_dtype_node, 3, label, disp) !mode=node, ndof=3
    end if
    if(iout(2) == 1) then
      label='velocity'
      call hecmw_result_add(hecmw_result_dtype_node, 3, label, vel) !mode=node, ndof=3
    end if
    if(iout(3) == 1) then
      label='acceleration'
      call hecmw_result_add(hecmw_result_dtype_node, 3, label, acc) !mode=node, ndof=3
    end if
    call hecmw_result_write_by_name(nameid)
    call hecmw_result_finalize()
    return
  end subroutine
 
  subroutine output_visfile(hecMESH, ifreq, disp, vel, acc, iout)
    !---- args
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    integer(kind=kint), intent(in)       :: ifreq
    real(kind=kreal), intent(in)         :: disp(:) !intend (numnodeDOF)
    real(kind=kreal), intent(in)         :: vel(:) !intend (numnodeDOF)
    real(kind=kreal), intent(in)         :: acc(:) !intend (numnodeDOF)
    integer(kind=kint), intent(in)       :: iout(3)
    !---- vals
    type(hecmwst_result_data)      :: fstrRESULT
    character(len=HECMW_NAME_LEN)  :: label
    integer(kind=kint)             :: ncomp, i
    !---- body
    ncomp = 0
    do i=1, 3
      if(iout(i) == 1) then
        ncomp = ncomp + 1
      end if
    end do
 
    call fstr_freq_result_init(hecmesh, ncomp, fstrresult)
    ncomp=1
    if(iout(1) == 1) then
      label = 'displace_abs'
      call fstr_freq_result_add(fstrresult, hecmesh, ncomp, 3, label, disp)
      ncomp = ncomp + 1
    end if
    if(iout(2) == 1) then
      label = 'velocity_abs'
      call fstr_freq_result_add(fstrresult, hecmesh, ncomp, 3, label, vel)
      ncomp = ncomp + 1
    end if
    if(iout(3) == 1) then
      label = 'acceleration_abs'
      call fstr_freq_result_add(fstrresult, hecmesh, ncomp, 3, label, acc)
      ncomp = ncomp + 1
    end if
 
    call fstr2hecmw_mesh_conv(hecmesh)
    call hecmw_visualize_init
    call hecmw_visualize( hecmesh, fstrresult, ifreq )
    call hecmw2fstr_mesh_conv(hecmesh)
    call hecmw_result_free(fstrresult)
  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
    integer(kind=kint)            :: NOD(NN)
    real(kind=kreal)              :: elecoord(3, nn), localcoord(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 )
      do i=1,nsur
        vect(3*nod(i)-2)=vect(3*nod(i)-2)+wg*h(i)*force(1)
        vect(3*nod(i)-1)=vect(3*nod(i)-1)+wg*h(i)*force(2)
        vect(3*nod(i)  )=vect(3*nod(i)  )+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
    !---- body
 
 
    fstrsolid%dunode = 0.d0
    call fstr_stiffmatrix( hecmesh, hecmat, fstrsolid, 0.d0, 0.d0 )
    call fstr_addbc(1, hecmesh, hecmat, fstrsolid, fstrparam, heclagmat, 2)
 
    call setmass(fstrsolid, hecmesh, hecmat, fstreig)
 
  end subroutine
 
  subroutine scaleeigenvector(fstrEIG, ntotaldof, nmode, eigenvector)
    !---- args
    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
    real(kind=kreal)   :: mas
    !---- body
 
    do imode=1,nmode
      mas = 0.0d0
      do idof=1,ntotaldof
        mas = mas + fstreig%mass(idof)*eigenvector(idof,imode)**2
      end do
      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(a, b, c)
    !---- args
    real(kind=kreal), intent(in)    :: a(:)
    real(kind=kreal), intent(in)    :: b(:)
    real(kind=kreal), intent(inout) :: c
    !---- vals
    !---- body
    c = dot_product(a, b)
    !Next, we need allreduce operation to implement distribute mesh.
    return
  end subroutine
 
  subroutine calcfreqcoeff(freqData, loadRe, loadIm, inpOmega, bjRe, bjIm)
    !---- args
    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(freqdata%eigVector(:,imode), loadre, ujfr)
      call calcdotproduct(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(cmplx(0.0d0, omega*time))
 
    do imode=1, freqdata%numMode
      s = size(freqdata%eigvector(:,imode))
      b = cmplx(bjre(imode), bjim(imode)) * a
      do idf=1, s
        c = b*cmplx(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 = cmplx(0.0d0, 1.0d0)*cmplx(omega, 0.0d0)*exp(cmplx(0.0d0, omega*time))
 
    do imode=1, freqdata%numMode
      s = size(freqdata%eigvector(:,imode))
      b = cmplx(bjre(imode), bjim(imode)) * a
      do idf=1, s
        c = b*cmplx(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 = cmplx(-1.0d0, 0.0d0)*cmplx(omega**2, 0.0d0)*exp(cmplx(0.0d0, omega*time))
 
    do imode=1, freqdata%numMode
      s = size(freqdata%eigvector(:,imode))
      b = cmplx(bjre(imode), bjim(imode)) * a
      do idf=1, s
        c = b*cmplx(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
 
  subroutine outputdyna_resfile(hecMESH, time, istp, dispre, dispim, velre, velim, accre, accim, iout)
    !---- args
    type(hecmwst_local_mesh), intent(in) :: hecMESH
    real(kind=kreal), intent(in)   :: time
    integer(kind=kint), intent(in) :: istp
    real(kind=kreal), intent(in)   :: dispre(:) !intend (numnodeDOF)
    real(kind=kreal), intent(in)   :: dispim(:) !intend (numnodeDOF)
    real(kind=kreal), intent(in)   :: velre(:) !intend (numnodeDOF)
    real(kind=kreal), intent(in)   :: velim(:) !intend (numnodeDOF)
    real(kind=kreal), intent(in)   :: accre(:) !intend (numnodeDOF)
    real(kind=kreal), intent(in)   :: accim(:) !intend (numnodeDOF)
    integer(kind=kint), intent(in) :: iout(3)
    !---- vals
    integer(kind=kint)              :: im, s
    character(len=HECMW_HEADER_LEN) :: header
    character(len=HECMW_MSG_LEN)    :: comment
    character(len=HECMW_NAME_LEN)   :: label, nameid
    real(kind=kreal), allocatable   :: absval(:)
    !---- body
 
    s = size(dispre)
    allocate(absval(s))
 
    nameid='fstrDYNA'
    header='*fstrresult'
    comment='frequency_result'
 
    call hecmw_result_init(hecmesh, istp, header, comment)
 
    label = "time"
    absval(1) = time
    call hecmw_result_add(hecmw_result_dtype_global, 1, label, absval)
 
    if(iout(1) == 1) then
      label='displacement_real'
      call hecmw_result_add(hecmw_result_dtype_node, 3, label, dispre) !mode=node, ndof=3
      label='displacement_imag'
      call hecmw_result_add(hecmw_result_dtype_node, 3, label, dispim)
      label='displacement_abs'
      absval(:) = abs(cmplx(dispre(:), dispim(:)))
      call hecmw_result_add(hecmw_result_dtype_node, 3, label, absval)
    end if
 
    if(iout(2) == 1) then
      label='velocity_real'
      call hecmw_result_add(hecmw_result_dtype_node, 3, label, velre) !mode=node, ndof=3
      label='velocity_imag'
      call hecmw_result_add(hecmw_result_dtype_node, 3, label, velim)
      label='velocity_abs'
      absval(:) = abs(cmplx(velre(:), velim(:)))
      call hecmw_result_add(hecmw_result_dtype_node, 3, label, absval)
    end if
 
    if(iout(3) == 1) then
      label='acceleration_real'
      call hecmw_result_add(hecmw_result_dtype_node, 3, label, accre) !mode=node, ndof=3
      label='acceleration_imag'
      call hecmw_result_add(hecmw_result_dtype_node, 3, label, accim)
      label='acceleration_abs'
      absval(:) = abs(cmplx(velre(:), velim(:)))
      call hecmw_result_add(hecmw_result_dtype_node, 3, label, absval)
    end if
 
    call hecmw_result_write_by_name(nameid)
    call hecmw_result_finalize()
 
    deallocate(absval)
    return
  end subroutine
 
  subroutine outputdyna_visfile(hecMESH, istp, dispre, dispim, velre, velim, accre, accim, iout)
    !---- args
    type(hecmwst_local_mesh), intent(inout) :: hecmesh
    integer(kind=kint), intent(in)          :: istp
    real(kind=kreal), intent(in)            :: dispre(:)
    real(kind=kreal), intent(in)            :: dispim(:)
    real(kind=kreal), intent(in)            :: velre(:)
    real(kind=kreal), intent(in)            :: velim(:)
    real(kind=kreal), intent(in)            :: accre(:)
    real(kind=kreal), intent(in)            :: accim(:)
    integer(kind=kint), intent(in)          :: iout(3)
    !---- vals
    type(hecmwst_result_data)     :: fstrRESULT
    character(len=HECMW_NAME_LEN) :: label
    integer(kind=kint)            :: s, ncomp, i
    real(kind=kreal), allocatable :: absval(:)
    !---- body
 
    s = size(dispre)
    allocate(absval(s))
 
    ncomp = 0
    do i=1, 3
      if(iout(i) == 1) then
        ncomp = ncomp + 3  !re, im, abs
      end if
    end do
 
    call fstr_freq_result_init(hecmesh, ncomp, fstrresult)  !disp, vel, acc
 
    ncomp = 1
 
    if(iout(1) == 1) then
      label = 'displace_real'
      call fstr_freq_result_add(fstrresult, hecmesh, ncomp, 3, label, dispre)
      ncomp = ncomp + 1
 
      label = 'displace_imag'
      call fstr_freq_result_add(fstrresult, hecmesh, ncomp, 3, label, dispim)
      ncomp = ncomp + 1
 
      label = 'displace_abs'
      absval(:) = abs(cmplx(dispre(:), dispim(:)))
      call fstr_freq_result_add(fstrresult, hecmesh, ncomp, 3, label, absval)
      ncomp = ncomp + 1
    end if
 
    if(iout(2) == 1) then
      label = 'velocity_real'
      call fstr_freq_result_add(fstrresult, hecmesh, ncomp, 3, label, velre)
      ncomp = ncomp + 1
 
      label = 'velocity_imag'
      call fstr_freq_result_add(fstrresult, hecmesh, ncomp, 3, label, velim)
      ncomp = ncomp + 1
 
      label = 'velocity_abs'
      absval(:) = abs(cmplx(velre(:), velim(:)))
      call fstr_freq_result_add(fstrresult, hecmesh, ncomp, 3, label, absval)
      ncomp = ncomp + 1
    end if
 
    if(iout(3) == 1) then
      label = 'acceleration_real'
      call fstr_freq_result_add(fstrresult, hecmesh, ncomp, 3, label, accre)
      ncomp = ncomp + 1
 
      label = 'acceleration_imag'
      call fstr_freq_result_add(fstrresult, hecmesh, ncomp, 3, label, accim)
      ncomp = ncomp + 1
 
      label = 'acceleration_abs'
      absval(:) = abs(cmplx(accre(:), accim(:)))
      call fstr_freq_result_add(fstrresult, hecmesh, ncomp, 3, label, absval)
      ncomp = ncomp + 1
    end if
 
    call fstr2hecmw_mesh_conv(hecmesh)
    call hecmw_visualize_init
    call hecmw_visualize( hecmesh, fstrresult, istp )
    call hecmw2fstr_mesh_conv(hecmesh)
    call hecmw_result_free(fstrresult)
 
    deallocate(absval)
  end subroutine
 
end module