fstr_EIG_output.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_eig_output
 
contains
 
  subroutine fstr_eigen_output(hecMESH, hecMAT, fstrEIG)
    use m_fstr
    use m_fstr_eig_lanczos_util
    use hecmw_solver_las
    implicit none
    type(hecmwst_local_mesh) :: hecMESH
    type(hecmwst_matrix)     :: hecMAT
    type(fstr_eigen)         :: fstrEIG
 
    integer(kind=kint) :: N, NP, NDOF, NNDOF, NPNDOF
    integer(kind=kint) :: i, j, k, in, jn, nget
    real(kind=kreal)   :: chk, gm
    real(kind=kreal), allocatable :: s(:), t(:), u(:), r(:)
    real(kind=kreal), pointer     :: mass(:), eigval(:), eigvec(:,:)
 
    n      = hecmat%N
    np     = hecmat%NP
    ndof   = hecmesh%n_dof
    nndof  = n *ndof
    npndof = np*ndof
 
    nget   =  fstreig%nget
    mass   => fstreig%mass
    eigval => fstreig%eigval
    eigvec => fstreig%eigvec
 
    allocate(fstreig%effmass(ndof*nget))
    allocate(fstreig%partfactor(ndof*nget))
    allocate(s(npndof))
    allocate(t(npndof))
    allocate(u(npndof))
    allocate(r(ndof))
    fstreig%effmass    = 0.0d0
    fstreig%partfactor = 0.0d0
 
    do i = 1, nget
      r  = 0.0d0
      gm = 0.0d0
      do j = 1, n
        do k = 1, ndof
          in = ndof*(j-1) + k
          r(k) = r(k) + mass(in)*eigvec(in,i)
          gm   = gm   + mass(in)*eigvec(in,i)*eigvec(in,i)
        enddo
      enddo
 
      call hecmw_allreduce_r(hecmesh, r, ndof, hecmw_sum)
      call hecmw_allreduce_r1(hecmesh, gm, hecmw_sum)
 
      gm = 1.0d0/gm
      do j = 1, ndof
        in = ndof*(i-1) + j
        fstreig%partfactor(in) = gm*r(j)
        fstreig%effmass(in)    = gm*r(j)*r(j)
      enddo
    enddo
 
    call eglist(hecmesh, hecmat, fstreig)
 
    if(myrank == 0)then
      write(imsg,*) ''
      write(imsg,*) '*----------------------------------------------*'
      write(imsg,*) '*--E I G E N V A L U E  C O N V E R G E N C E--*'
      write(imsg,*) 'Absolute residual =  |(||Kx - lambda*Mx||)|'
      write(imsg,*) '   Iter.#   Eigenvalue    Abs. Residual  '
      write(ilog,*) '   Iter.#   Eigenvalue    Abs. Residual  '
      write(imsg,*) '   *-----*  *---------*  *--------------*'
    endif
 
    do j = 1, fstreig%nget
      do i = 1, nndof
        u(i) = eigvec(i,j)
      enddo
      call hecmw_matvec(hecmesh, hecmat, u, t)
 
      s = 0.0d0
      do i = 1, nndof
        s(i) = mass(i)*eigvec(i,j)
      enddo
 
      chk = 0.0d0
      do i = 1,nndof
        chk = chk + (t(i) - (eigval(j)-fstreig%sigma)*s(i))**2
      enddo
      call hecmw_allreduce_r1(hecmesh, chk, hecmw_sum)
      chk = dsqrt(chk)
 
      if(myrank == 0)then
        write(imsg,'(2x,i5,2x,1p5e15.6)') j, eigval(j), chk
        write(ilog,'(i5,1p5e12.4)') j, eigval(j), chk
      endif
    enddo
 
    if(myrank == 0)then
      write(imsg,*)'*    ---END Eigenvalue listing---     *'
    endif
 
    deallocate(s)
    deallocate(t)
    deallocate(u)
    deallocate(r)
  end subroutine fstr_eigen_output
 
  subroutine fstr_eigen_make_result(hecMESH, hecMAT, fstrEIG, fstrRESULT)
    use m_fstr
    use m_hecmw2fstr_mesh_conv
    use hecmw_util
    implicit none
    type(hecmwst_local_mesh)  :: hecmesh
    type(hecmwst_matrix)      :: hecMAT
    type(fstr_eigen)          :: fstrEIG
    type(hecmwst_result_data) :: fstrRESULT
 
    integer(kind=kint) :: i, istep, nget, NP, NDOF, NPNDOF, totalmpc, MPC_METHOD
    real(kind=kreal)   :: t1
    real(kind=kreal), pointer :: eigvec(:,:)
    real(kind=kreal), allocatable :: x(:), egval(:)
    real(kind=kreal), allocatable :: disp3(:), rot3(:)
    character(len=HECMW_HEADER_LEN) :: header
    character(len=HECMW_MSG_LEN)    :: comment
    character(len=HECMW_NAME_LEN)   :: label
    character(len=HECMW_NAME_LEN)   :: nameID
 
    nget   = fstreig%nget
    np     = hecmat%NP
    ndof   = hecmat%NDOF
    npndof = hecmat%NP*hecmat%NDOF
    !totalmpc = hecMESH%mpc%n_mpc
    !call hecmw_allreduce_I1 (hecMESH, totalmpc, hecmw_sum)
 
    eigvec => fstreig%eigvec
 
    allocate(x(npndof))
    x = 0.0d0
    allocate(egval(1))
 
    do istep = 1, nget
      egval(1) = fstreig%eigval(istep)
      do i=1,npndof
        x(i) = eigvec(i,istep)
      enddo
 
      !if (totalmpc > 0) then
      !  MPC_METHOD = hecmw_mat_get_mpc_method(hecMAT)
      !  if (MPC_METHOD < 1 .or. 3 < MPC_METHOD) MPC_METHOD = 3
      !  if (MPC_METHOD == 3) then  ! elimination
      !    call hecmw_tback_x_33(hecMESH, X, t1)
      !  else
      !    if (hecMESH%my_rank.eq.0) write(0,*) "### ERROR: MPC_METHOD must set to 3"
      !    stop
      !  endif
      !endif
 
      call hecmw_update_r(hecmesh, x, hecmat%NP, ndof)
 
      ! For 6-DOF elements (e.g. shells 731/741), the eigenvector holds 3
      ! translational + 3 rotational components per node. Writing all 6 as
      ! DISPLACEMENT makes VTK/ParaView interpret it as a symmetric tensor
      ! (NumberOfComponents=6), breaking Warp By Vector. Mirror the static
      ! solver (make_result.f90, ndof==6) and split DISPLACEMENT (1:3) from
      ! ROTATION (4:6). Other NDOF (e.g. solid 3) keep the pass-through.
      if( ndof == 6 ) then
        allocate(disp3(np*3))
        allocate(rot3(np*3))
        disp3 = 0.0d0
        rot3  = 0.0d0
        do i = 1, np
          disp3((i-1)*3+1:(i-1)*3+3) = x((i-1)*ndof+1:(i-1)*ndof+3)
          rot3((i-1)*3+1:(i-1)*3+3) = x((i-1)*ndof+4:(i-1)*ndof+6)
        enddo
      endif
 
      if( iresult.eq.1 ) then
        header = "*fstrresult"
        comment = "eigen_result"
        call hecmw_result_init(hecmesh,istep,header,comment)
        label = "EIGENVALUE"
        call hecmw_result_add(hecmw_result_dtype_global,1,label,egval)
        if( ndof == 6 ) then
          label = "DISPLACEMENT"
          call hecmw_result_add(hecmw_result_dtype_node,3,label,disp3)
          label = "ROTATION"
          call hecmw_result_add(hecmw_result_dtype_node,3,label,rot3)
        else
          label = "DISPLACEMENT"
          call hecmw_result_add(hecmw_result_dtype_node,ndof,label,x)
        endif
        nameid = "fstrRES"
        call hecmw_result_write_by_name(nameid)
        call hecmw_result_finalize
      endif
 
      if( ivisual.eq.1 ) then
        call hecmw_nullify_result_data(fstrresult)
        fstrresult%ng_component = 1
        fstrresult%ne_component = 0
        allocate(fstrresult%ng_dof(1))
        allocate(fstrresult%global_label(1))
        allocate(fstrresult%global_val_item(1))
        fstrresult%ng_dof(1) = 1
        fstrresult%global_label(1) = 'EIGENVALUE'
        fstrresult%global_val_item(1) = egval(1)
        if( ndof == 6 ) then
          fstrresult%nn_component = 2
          allocate(fstrresult%nn_dof(2))
          allocate(fstrresult%node_label(2))
          allocate(fstrresult%node_val_item(6*np))
          fstrresult%nn_dof(1) = 3
          fstrresult%nn_dof(2) = 3
          fstrresult%node_label(1) = 'DISPLACEMENT'
          fstrresult%node_label(2) = 'ROTATION'
          do i = 1, np
            fstrresult%node_val_item((i-1)*6+1:(i-1)*6+3) = disp3((i-1)*3+1:(i-1)*3+3)
            fstrresult%node_val_item((i-1)*6+4:(i-1)*6+6) = rot3((i-1)*3+1:(i-1)*3+3)
          enddo
        else
          fstrresult%nn_component = 1
          allocate(fstrresult%nn_dof(1))
          allocate(fstrresult%node_label(1))
          allocate(fstrresult%node_val_item(ndof*np))
          fstrresult%nn_dof(1) = ndof
          fstrresult%node_label(1) = 'DISPLACEMENT'
          fstrresult%node_val_item = x
        endif
        call fstr2hecmw_mesh_conv(hecmesh)
        call hecmw_visualize_init
        call hecmw_visualize( hecmesh, fstrresult, istep )
        call hecmw_visualize_finalize
        call hecmw2fstr_mesh_conv(hecmesh)
        call hecmw_result_free(fstrresult)
      endif
 
      if( ndof == 6 ) then
        deallocate(disp3)
        deallocate(rot3)
      endif
    enddo
 
    deallocate(x)
 
  end subroutine fstr_eigen_make_result
 
  subroutine eglist(hecMESH, hecMAT, fstrEIG)
    use m_fstr
    use hecmw_util
    use m_fstr_eig_lanczos_util
    implicit none
    type(hecmwst_local_mesh) :: hecMESH
    type(hecmwst_matrix)     :: hecMAT
    type(fstr_eigen)         :: fstrEIG
 
    integer(kind=kint) :: NDOF
    integer(kind=kint) :: i, j, in, iter ,nget
    real(kind=kreal)   :: pi, angle, freq, pf(3), em(3)
    real(kind=kreal), pointer :: eigval(:)
 
    ndof = hecmat%NDOF
    nget = fstreig%nget
    iter = fstreig%iter
    pi   = 4.0d0 * datan(1.0d0)
    eigval => fstreig%eigval
 
    if(myrank == 0)then
      write(ilog,*)""
      write(ilog,"(a)")"********************************"
      write(ilog,"(a)")"*RESULT OF EIGEN VALUE ANALYSIS*"
      write(ilog,"(a)")"********************************"
      write(ilog,"(a)")""
      write(ilog,"(a,i8)")"NUMBER OF ITERATIONS = ",iter
      write(ilog,"(a,1pe12.4)")"TOTAL MASS = ",fstreig%totalmass
      write(ilog,"(a)")""
      write(ilog,"(3a)")"                   ANGLE       FREQUENCY   ",&
        "PARTICIPATION FACTOR                EFFECTIVE MASS"
      write(ilog,"(3a)")"  NO.  EIGENVALUE  FREQUENCY   (HZ)        ",&
        "X           Y           Z           X           Y           Z"
      write(ilog,"(3a)")"  ---  ----------  ----------  ----------  ",&
        "----------  ----------  ----------  ----------  ----------  ----------"
      write(*,*)""
      write(*,"(a)")"#----------------------------------#"
      write(*,"(a)")"#  RESULT OF EIGEN VALUE ANALYSIS  #"
      write(*,"(a)")"#----------------------------------#"
      write(*,"(a)")""
      write(*,"(a,i8)")"### NUMBER OF ITERATIONS = ",iter
      write(*,"(a,1pe12.4)")"### TOTAL MASS = ",fstreig%totalmass
      write(*,"(a)")""
      write(*,"(3a)")"       PERIOD     FREQUENCY  ",&
        "PARTICIPATION FACTOR             EFFECTIVE MASS"
      write(*,"(3a)")"  NO.  [Sec]      [HZ]       ",&
        "X          Y          Z          X          Y          Z"
      write(*,"(3a)")"  ---  ---------  ---------  ",&
        "---------  ---------  ---------  ---------  ---------  ---------"
 
      in = 0
      do i = 1, nget
        in = in + 1
        if(eigval(i) < 0.0d0) eigval(i) = 0.0d0
        angle = dsqrt(eigval(i))
        freq = angle*0.5d0/pi
 
        pf = 0.0d0
        em = 0.0d0
        do j = 1, min(ndof, 3)
          pf(j) = fstreig%partfactor(ndof*(i-1) + j)
          em(j) = fstreig%effmass(ndof*(i-1) + j)
        enddo
 
        write(ilog,'(I5,1P9E12.4)') in, eigval(i), angle, freq, pf(1), pf(2), pf(3), em(1), em(2), em(3)
        write(*   ,'(I5,1P8E11.3)') in, 1.0d0/freq, freq, pf(1), pf(2), pf(3), em(1), em(2), em(3)
      enddo
      write(ilog,*)""
      write(*,*)""
    endif
 
  end subroutine eglist
end module m_fstr_eig_output