fstr_Residual.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_residual
  use hecmw
  use m_fstr_elemact
  implicit none
 
  public :: fstr_update_ndforce
  public :: fstr_update_ndforce_spc
  public :: fstr_update_reaction_spc
  public :: fstr_get_residual
  public :: fstr_get_norm_contact
  public :: fstr_assemble_residual_contact
  public :: fstr_get_norm_para_contact
  public :: fstr_get_x_norm_contact
  public :: fstr_get_potential
 
  private :: fstr_update_ndforce_mpc
 
contains
 
  !C---------------------------------------------------------------------*
  subroutine fstr_update_ndforce(cstep,hecMESH,hecMAT,fstrSOLID,conMAT)
    !C---------------------------------------------------------------------*
    use m_fstr
    use muload
    use m_fstr_spring
    integer(kind=kint), intent(in)       :: cstep
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    type(hecmwst_matrix), intent(inout)  :: hecmat
    type(fstr_solid), intent(inout)      :: fstrsolid
    type(hecmwst_matrix), intent(inout), optional  :: conmat
    !    Local variables
    integer(kind=kint) :: ndof, idof
    real(kind=kreal)   :: factor
 
    factor = fstrsolid%factor(2)
 
    !    Set residual load
    do idof=1, hecmesh%n_node*  hecmesh%n_dof
      hecmat%B(idof)=fstrsolid%GL(idof)-fstrsolid%QFORCE(idof)
    end do
    ndof = hecmat%NDOF
 
    call fstr_update_ndforce_spring( cstep, hecmesh, fstrsolid, hecmat%B )
 
    !    Consider Uload
    call uresidual( cstep, factor, hecmat%B )
 
    !    Consider EQUATION condition
    call fstr_update_ndforce_mpc( hecmesh, hecmat%B )
 
    !    Consider SPC condition
    call fstr_update_ndforce_spc( cstep, hecmesh, fstrsolid, hecmat%B )
    if(present(conmat)) call fstr_update_ndforce_spc( cstep, hecmesh, fstrsolid, conmat%B )
 
    !
    call hecmw_update_r(hecmesh,hecmat%B,hecmesh%n_node, ndof)
  end subroutine fstr_update_ndforce
 
  subroutine fstr_update_ndforce_mpc( hecMESH, B )
    use m_fstr
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    real(kind=kreal), intent(inout)      :: b(:) 
    !    Local variables
    integer(kind=kint) ndof, ig0, is0, ie0, ik, in, idof
    real(kind=kreal) :: rhs, lambda
 
    ndof = hecmesh%n_dof
    outer: do ig0=1,hecmesh%mpc%n_mpc
      is0= hecmesh%mpc%mpc_index(ig0-1)+1
      ie0= hecmesh%mpc%mpc_index(ig0)
      do ik= is0, ie0
        if (hecmesh%mpc%mpc_dof(ik) > ndof) cycle outer
      enddo
      ! Suppose the lagrange multiplier= first dof of first node
      in = hecmesh%mpc%mpc_item(is0)
      idof = hecmesh%mpc%mpc_dof(is0)
      rhs = hecmesh%mpc%mpc_val(is0)
      lambda = b(ndof*(in-1)+idof)/rhs
      ! update nodal residual
      do ik= is0, ie0
        in = hecmesh%mpc%mpc_item(ik)
        idof = hecmesh%mpc%mpc_dof(ik)
        rhs = hecmesh%mpc%mpc_val(ik)
        b(ndof*(in-1)+idof) = b(ndof*(in-1)+idof) - rhs*lambda
      enddo
    enddo outer
  end subroutine fstr_update_ndforce_mpc
 
  subroutine fstr_update_ndforce_spc( cstep, hecMESH, fstrSOLID, B )
    use m_fstr
    integer(kind=kint), intent(in)       :: cstep
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    type(fstr_solid), intent(in)         :: fstrsolid
    real(kind=kreal), intent(inout)      :: b(:) 
    !    Local variables
    integer(kind=kint) ndof, ig0, ig, ityp, is0, ie0, ik, in, idof1, idof2, idof
    integer(kind=kint) :: grpid
 
    ndof = hecmesh%n_dof
 
    !    Clear RHS at constrained DOFs (SPC)
    do ig0= 1, fstrsolid%BOUNDARY_ngrp_tot
      grpid = fstrsolid%BOUNDARY_ngrp_GRPID(ig0)
      if( .not. fstr_isboundaryactive( fstrsolid, grpid, cstep ) ) cycle
      ig= fstrsolid%BOUNDARY_ngrp_ID(ig0)
      ityp= fstrsolid%BOUNDARY_ngrp_type(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)
        idof1 = ityp/10
        idof2 = ityp - idof1*10
        if( fstrsolid%BOUNDARY_ngrp_rotID(ig0) > 0 ) then
          idof1 = 1
          idof2 = ndof
        end if
        do idof=idof1,idof2
          b( ndof*(in-1) + idof ) = 0.d0
        enddo
      enddo
    enddo
  end subroutine fstr_update_ndforce_spc
 
  subroutine fstr_update_reaction_spc( cstep, hecMESH, fstrSOLID )
    use m_fstr
    integer(kind=kint), intent(in)       :: cstep
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    type(fstr_solid), intent(in)         :: fstrsolid
    !    Local variables
    integer(kind=kint) :: ndof, ig0, ig, ityp, is0, ie0, ik, in, idof1, idof2, idof
    integer(kind=kint) :: grpid
 
    ndof = hecmesh%n_dof
    fstrsolid%REACTION = 0.d0
 
    ! displacement boundary
    do ig0= 1, fstrsolid%BOUNDARY_ngrp_tot
      grpid = fstrsolid%BOUNDARY_ngrp_GRPID(ig0)
      if( .not. fstr_isboundaryactive( fstrsolid, grpid, cstep ) ) cycle
      ig= fstrsolid%BOUNDARY_ngrp_ID(ig0)
      ityp= fstrsolid%BOUNDARY_ngrp_type(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)
        idof1 = ityp/10
        idof2 = ityp - idof1*10
        if( fstrsolid%BOUNDARY_ngrp_rotID(ig0) > 0 ) then
          idof1 = 1
          idof2 = ndof
        end if
        do idof=idof1,idof2
          fstrsolid%REACTION(ndof*(in-1)+idof) = fstrsolid%QFORCE(ndof*(in-1)+idof)
          !count embed force as reaction force
          if( associated(fstrsolid%EMBED_NFORCE) ) fstrsolid%REACTION(ndof*(in-1)+idof) = &
                &  fstrsolid%QFORCE(ndof*(in-1)+idof) - fstrsolid%EMBED_NFORCE(ndof*(in-1)+idof)
        enddo
      enddo
    enddo
 
    ! velocity boundary (dynamic only; following type constrains the DOF)
    if( fstrsolid%VELOCITY_type /= kbcinitial ) then
      do ig0= 1, fstrsolid%VELOCITY_ngrp_tot
        grpid = fstrsolid%VELOCITY_ngrp_GRPID(ig0)
        if( .not. fstr_isboundaryactive( fstrsolid, grpid, cstep ) ) cycle
        ig= fstrsolid%VELOCITY_ngrp_ID(ig0)
        ityp= fstrsolid%VELOCITY_ngrp_type(ig0)
        is0= hecmesh%node_group%grp_index(ig-1) + 1
        ie0= hecmesh%node_group%grp_index(ig  )
        idof1 = ityp/10
        idof2 = ityp - idof1*10
        do ik= is0, ie0
          in   = hecmesh%node_group%grp_item(ik)
          do idof=idof1,idof2
            fstrsolid%REACTION(ndof*(in-1)+idof) = fstrsolid%QFORCE(ndof*(in-1)+idof)
          enddo
        enddo
      enddo
    endif
 
    ! acceleration boundary (dynamic only; following type constrains the DOF)
    if( fstrsolid%ACCELERATION_type /= kbcinitial ) then
      do ig0= 1, fstrsolid%ACCELERATION_ngrp_tot
        grpid = fstrsolid%ACCELERATION_ngrp_GRPID(ig0)
        if( .not. fstr_isboundaryactive( fstrsolid, grpid, cstep ) ) cycle
        ig= fstrsolid%ACCELERATION_ngrp_ID(ig0)
        ityp= fstrsolid%ACCELERATION_ngrp_type(ig0)
        is0= hecmesh%node_group%grp_index(ig-1) + 1
        ie0= hecmesh%node_group%grp_index(ig  )
        idof1 = ityp/10
        idof2 = ityp - idof1*10
        do ik= is0, ie0
          in   = hecmesh%node_group%grp_item(ik)
          do idof=idof1,idof2
            fstrsolid%REACTION(ndof*(in-1)+idof) = fstrsolid%QFORCE(ndof*(in-1)+idof)
          enddo
        enddo
      enddo
    endif
  end subroutine fstr_update_reaction_spc
 
 
  subroutine fstr_calc_residual_vector(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM)
    use m_fstr_update
    use m_fstr_ass_load
    implicit none
    integer, intent(in)                   :: cstep
    type (hecmwST_local_mesh)             :: hecMESH
    type (hecmwST_matrix)                 :: hecMAT
    type (fstr_solid)                     :: fstrSOLID
    real(kind=kreal), intent(in)          :: ctime     
    real(kind=kreal), intent(in)          :: dtime     
    type (fstr_param)                     :: fstrPARAM
    real(kind=kreal), intent(in) :: tincr
    integer(kind=kint) :: iter
 
    ! ----- update the strain, stress, and internal force
    call fstr_updatenewton(hecmesh, hecmat, fstrsolid, ctime, tincr, iter)
 
    ! ----- Set residual
    if( fstrsolid%DLOAD_follow /= 0 .or. fstrsolid%CLOAD_ngrp_rot /= 0 ) &
      & call fstr_ass_load(cstep, ctime+dtime, hecmesh, hecmat, fstrsolid, fstrparam )
 
 
    if( fstrsolid%elemact%ELEMACT_egrp_tot > 0 ) then
      call fstr_update_elemact_solid_by_value( hecmesh, fstrsolid, cstep, ctime )
    endif
    
    call fstr_update_ndforce(cstep, hecmesh, hecmat, fstrsolid)
 
  end subroutine fstr_calc_residual_vector
 
  subroutine fstr_calc_residual_vector_with_x(hecMESH, hecMAT, fstrSOLID, ctime, tincr, iter, cstep, dtime, fstrPARAM)
    use m_fstr_update, only: fstr_updatenewton
    use m_fstr_nodalkinematics, only: fstr_apply_solution_increment
    implicit none
    integer, intent(in)                   :: cstep
    type (hecmwST_local_mesh)             :: hecMESH
    type (hecmwST_matrix)                 :: hecMAT
    type (fstr_solid)                     :: fstrSOLID
    real(kind=kreal), intent(in)          :: ctime     
    real(kind=kreal), intent(in)          :: dtime     
    type (fstr_param)                     :: fstrPARAM
    real(kind=kreal), intent(in) :: tincr
    integer(kind=kint) :: iter
 
    integer :: i
    real(kind=kreal) :: dunode_bak(hecmat%ndof*hecmesh%n_node)
    real(kind=kreal), allocatable :: shell_dtriad_bak(:), shell_ddrill_bak(:)
 
    do i=1, hecmat%ndof*hecmesh%n_node
      dunode_bak(i) = fstrsolid%dunode(i)
    enddo
    if( associated(fstrsolid%shell_dtriad) ) then
      allocate(shell_dtriad_bak(size(fstrsolid%shell_dtriad)))
      shell_dtriad_bak(:) = fstrsolid%shell_dtriad(:)
    endif
    if( associated(fstrsolid%shell_ddrill) ) then
      allocate(shell_ddrill_bak(size(fstrsolid%shell_ddrill)))
      shell_ddrill_bak(:) = fstrsolid%shell_ddrill(:)
    endif
    call fstr_apply_solution_increment( hecmesh, fstrsolid, hecmat%ndof, hecmat%X )
    call fstr_calc_residual_vector(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam)
    do i=1, hecmat%ndof*hecmesh%n_node
      fstrsolid%dunode(i) = dunode_bak(i)
    enddo
    if( allocated(shell_dtriad_bak) ) then
      fstrsolid%shell_dtriad(:) = shell_dtriad_bak(:)
      deallocate(shell_dtriad_bak)
    endif
    if( allocated(shell_ddrill_bak) ) then
      fstrsolid%shell_ddrill(:) = shell_ddrill_bak(:)
      deallocate(shell_ddrill_bak)
    endif
  end subroutine fstr_calc_residual_vector_with_x
 
  real(kind=kreal) function fstr_get_residual( force, hecMESH )
    use m_fstr
    real(kind=kreal), intent(in)         :: force(:)
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    integer :: ndof
    ndof = hecmesh%n_dof
    call hecmw_innerproduct_r(hecmesh,ndof,force,force,fstr_get_residual)
  end function
 
  real(kind=kreal) function fstr_get_energy( force, displacement, hecMESH )
    use m_fstr
    real(kind=kreal), intent(in)         :: force(:), displacement(:)
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    integer :: ndof
    ndof = hecmesh%n_dof
    call hecmw_innerproduct_r(hecmesh, ndof, force, displacement, fstr_get_energy)
  end function
 
  real(kind=kreal) function fstr_get_norm_contact(flag,hecMESH,hecMAT,fstrSOLID,hecLagMAT)
    use m_fstr
    type(hecmwst_local_mesh), intent(in)             :: hecmesh
    type(hecmwst_matrix), intent(in)                 :: hecmat
    type(fstr_solid), intent(in)                     :: fstrsolid
    type(hecmwst_matrix_lagrange), intent(in)        :: heclagmat
    character(len=13)                                :: flag
    real(kind=kreal) :: tmp1, tmp2, bi
    integer :: i, i0, ndof
    if( flag=='residualForce' )then
      ndof = hecmesh%n_dof
      call hecmw_innerproduct_r(hecmesh,ndof,hecmat%B,hecmat%B,tmp1)
      tmp2 = 0.0d0
      i0 = hecmesh%n_node*ndof
      do i=1,heclagmat%num_lagrange
        bi = hecmat%B(i0+i)
        tmp2 = tmp2 + bi*bi
      enddo
      call hecmw_allreduce_r1(hecmesh,tmp2,hecmw_sum)
      fstr_get_norm_contact = tmp1 + tmp2
    elseif( flag=='        force' )then
      call hecmw_innerproduct_r(hecmesh,ndof,fstrsolid%QFORCE,fstrsolid%QFORCE,fstr_get_norm_contact)
    endif
  end function
 
  subroutine fstr_assemble_residual_contact(hecMAT, hecLagMAT, conMAT, hecMESH, resid_vec, nresid)
    use m_fstr
    implicit none
    type(hecmwst_matrix), intent(in)          :: hecmat
    type(hecmwst_matrix_lagrange), intent(in) :: heclagmat
    type(hecmwst_matrix), intent(in)          :: conmat
    type(hecmwst_local_mesh), intent(in)      :: hecmesh
    real(kind=kreal), intent(out)             :: resid_vec(:)
    integer(kind=kint), intent(out)           :: nresid
 
    integer(kind=kint) :: i, ndof, nndof, npndof, num_lagrange
 
    ndof = conmat%ndof
    nndof = conmat%N * ndof
    npndof = conmat%NP * ndof
    num_lagrange = heclagmat%num_lagrange
    nresid = nndof + num_lagrange
 
    ! Copy conMAT%B and assemble across MPI
    do i = 1, npndof
      resid_vec(i) = conmat%B(i)
    enddo
    call hecmw_assemble_r(hecmesh, resid_vec, conmat%NP, conmat%NDOF)
 
    ! Add hecMAT%B (node DOFs only)
    do i = 1, nndof
      resid_vec(i) = resid_vec(i) + hecmat%B(i)
    enddo
 
    ! Copy Lagrange multiplier residual
    do i = 1, num_lagrange
      resid_vec(nndof + i) = conmat%B(npndof + i)
    enddo
 
  end subroutine fstr_assemble_residual_contact
 
  !
  function fstr_get_norm_para_contact(hecMAT,hecLagMAT,conMAT,hecMESH) result(rhsB)
    use m_fstr
    implicit none
    type(hecmwst_matrix), intent(in)                 :: hecmat
    type(hecmwst_matrix_lagrange), intent(in)        :: heclagmat
    type(hecmwst_matrix), intent(in)                 :: conmat
    type(hecmwst_local_mesh), intent(in)             :: hecmesh
    !
    real(kind=kreal) ::  rhsb
    integer(kind=kint) ::  i,ndof,nndof,npndof,num_lagrange
    real(kind=kreal), allocatable   :: rhs_con(:)
    real(kind=kreal), pointer :: rhs_lag(:)
 
    ndof = conmat%ndof
    nndof = conmat%N * ndof
    npndof = conmat%NP * ndof
    num_lagrange = heclagmat%num_lagrange
 
    allocate(rhs_con(npndof))
    do i=1,npndof
      rhs_con(i) = conmat%B(i)
    enddo
    call hecmw_assemble_r(hecmesh, rhs_con, conmat%NP, conmat%NDOF)
 
    do i=1,nndof
      rhs_con(i) = rhs_con(i) + hecmat%B(i)
    enddo
 
    rhs_lag => conmat%B(npndof+1:npndof+num_lagrange)
 
    rhsb = dot_product(rhs_con(1:nndof), rhs_con(1:nndof)) + dot_product(rhs_lag(:), rhs_lag(:))
    call hecmw_allreduce_r1(hecmesh, rhsb, hecmw_sum)
    deallocate(rhs_con)
 
  end function fstr_get_norm_para_contact
 
  function fstr_get_x_norm_contact(hecMAT,hecLagMAT,hecMESH) result(rhsX)
    use m_fstr
    implicit none
    type(hecmwst_matrix), intent(in)                 :: hecmat
    type(hecmwst_matrix_lagrange), intent(in)        :: heclagmat
    type(hecmwst_local_mesh), intent(in)             :: hecmesh
    real(kind=kreal)   ::  rhsx
    integer(kind=kint) :: nndof, npndof, i
 
    nndof = hecmat%N * hecmat%NDOF
    npndof = hecmat%NP * hecmat%NDOF
    rhsx = 0.d0
    do i=1,nndof
      rhsx = rhsx + hecmat%X(i) ** 2
    end do
    do i=1,heclagmat%num_lagrange
      rhsx = rhsx + hecmat%X(npndof+i) ** 2
    end do
    call hecmw_allreduce_r1(hecmesh, rhsx, hecmw_sum)
 
  end function fstr_get_x_norm_contact
 
  !C---------------------------------------------------------------------*
  function fstr_get_potential(cstep,hecMESH,hecMAT,fstrSOLID,ptype) result(potential)
    use m_fstr
    use muload
    use m_fstr_spring
    integer(kind=kint), intent(in)       :: cstep
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    type(hecmwst_matrix), intent(inout)  :: hecmat
    type(fstr_solid), intent(inout)      :: fstrsolid
    integer(kind=kint), intent(in)       :: ptype
    real(kind=kreal)   ::  potential
    !    Local variables
    integer(kind=kint) :: ndof, i, icel, icel0, ig
    real(kind=kreal), allocatable :: totdisp(:), totload(:)
    real(kind=kreal) :: factor
    real(kind=kreal) :: extenal_work, internal_work
 
    ndof = hecmat%NDOF
    factor = fstrsolid%factor(2)
    potential = 0.d0
 
    allocate(totload(ndof*hecmesh%n_node))
 
    ! Set external load
    do i=1,ndof*hecmesh%n_node
      totload(i) = 0.5d0*(fstrsolid%GL(i)+fstrsolid%GL0(i))
    end do
    !    Consiter Spring
    call fstr_update_ndforce_spring( cstep, hecmesh, fstrsolid, totload )
    !    Consider Uload
    call uresidual( cstep, factor, totload )
    !    Consider EQUATION condition
    call fstr_update_ndforce_mpc( hecmesh, totload )
 
    if( ptype == 1 ) then 
      ! calc internal work
      internal_work = 0.d0
      do icel0=1,hecmesh%ne_internal
        icel=hecmesh%elem_internal_list(icel0)
        do ig=1,size(fstrsolid%elements(icel)%gausses)
          internal_work = internal_work + fstrsolid%elements(icel)%gausses(ig)%strain_energy
        enddo
      enddo
      call hecmw_allreduce_r1(hecmesh, internal_work, hecmw_sum)
      ! calc external work
      call hecmw_innerproduct_r(hecmesh,ndof,totload,fstrsolid%dunode,extenal_work)      
      potential = internal_work - extenal_work
    else if( ptype == 2 ) then ! multiply internal force with displacement
      ! calc internal work
      call hecmw_innerproduct_r(hecmesh,ndof,0.5d0*(fstrsolid%QFORCE+fstrsolid%QFORCE_bak),fstrsolid%dunode,internal_work)
      ! calc external work
      call hecmw_innerproduct_r(hecmesh,ndof,totload,fstrsolid%dunode,extenal_work)
      potential = internal_work - extenal_work
    else if( ptype == 3 ) then ! norm of residual
      ! Consider internal force
      totload(:) = fstrsolid%QFORCE(:) - totload(:)
      !    Consider SPC condition
      call fstr_update_ndforce_spc( cstep, hecmesh, fstrsolid, totload )
      ! calc norm of residual
      call hecmw_innerproduct_r(hecmesh,ndof,totload,totload,potential)
      potential = dsqrt(potential)
    else
      stop "ptype error in fstr_get_potential"
    endif
    !write(IMSG,*) "internal_work,extenal_work",internal_work,extenal_work
  end function fstr_get_potential
 
  function fstr_get_potential_with_x(cstep,hecMESH,hecMAT,fstrSOLID,ptype) result(potential)
    use m_fstr
    use muload
    use m_fstr_spring
    use m_fstr_nodalkinematics, only: fstr_apply_solution_increment
    implicit none
    integer(kind=kint), intent(in)       :: cstep
    type(hecmwst_local_mesh), intent(in) :: hecmesh
    type(hecmwst_matrix), intent(inout)  :: hecmat
    type(fstr_solid), intent(inout)      :: fstrsolid
    integer(kind=kint), intent(in)       :: ptype
    real(kind=kreal)   ::  potential
 
    integer :: i
    real(kind=kreal) :: dunode_bak(hecmat%ndof*hecmesh%n_node)
    real(kind=kreal), allocatable :: shell_dtriad_bak(:), shell_ddrill_bak(:)
 
    do i=1, hecmat%ndof*hecmesh%n_node
      dunode_bak(i) = fstrsolid%dunode(i)
    enddo
    if( associated(fstrsolid%shell_dtriad) ) then
      allocate(shell_dtriad_bak(size(fstrsolid%shell_dtriad)))
      shell_dtriad_bak(:) = fstrsolid%shell_dtriad(:)
    endif
    if( associated(fstrsolid%shell_ddrill) ) then
      allocate(shell_ddrill_bak(size(fstrsolid%shell_ddrill)))
      shell_ddrill_bak(:) = fstrsolid%shell_ddrill(:)
    endif
    call fstr_apply_solution_increment( hecmesh, fstrsolid, hecmat%ndof, hecmat%X )
    potential = fstr_get_potential(cstep,hecmesh,hecmat,fstrsolid,ptype)
    do i=1, hecmat%ndof*hecmesh%n_node
      fstrsolid%dunode(i) = dunode_bak(i)
    enddo
    if( allocated(shell_dtriad_bak) ) then
      fstrsolid%shell_dtriad(:) = shell_dtriad_bak(:)
      deallocate(shell_dtriad_bak)
    endif
    if( allocated(shell_ddrill_bak) ) then
      fstrsolid%shell_ddrill(:) = shell_ddrill_bak(:)
      deallocate(shell_ddrill_bak)
    endif
  end function fstr_get_potential_with_x
 
  subroutine plot_potential_graph( cstep, hecMESH, hecMAT, fstrSOLID, fstrPARAM, &
    restrt_step_num, sub_step, ctime, dtime, iter, tincr )
    use m_fstr
    use m_fstr_nodalkinematics, only: fstr_apply_solution_increment
    integer, intent(in)                   :: cstep
    type (hecmwST_local_mesh)             :: hecMESH
    type (hecmwST_matrix)                 :: hecMAT
    type (fstr_solid)                     :: fstrSOLID
    integer, intent(in)                   :: sub_step
    real(kind=kreal), intent(in)          :: ctime     
    real(kind=kreal), intent(in)          :: dtime     
    type (fstr_param)                     :: fstrPARAM
    type (hecmwST_matrix_lagrange)        :: hecLagMAT
    integer(kind=kint) :: restrt_step_num
 
    integer(kind=kint) :: iter, i, j
    integer(kind=kint), parameter :: ntot = 30
    real(kind=kreal)   :: tincr, alpha, dulen
    real(kind=kreal) :: pot(3)
    real(kind=kreal), allocatable :: dunode_bak(:)
    real(kind=kreal), allocatable :: shell_dtriad_bak(:), shell_ddrill_bak(:)
 
    dulen = dsqrt(dot_product(fstrsolid%dunode(:),fstrsolid%dunode(:)))
    write(imsg,*) "dulen",dulen
 
    allocate(dunode_bak(size(fstrsolid%dunode)))
    do i=1, hecmat%ndof*hecmesh%n_node
      dunode_bak(i) = fstrsolid%dunode(i)
    enddo
    if( associated(fstrsolid%shell_dtriad) ) then
      allocate(shell_dtriad_bak(size(fstrsolid%shell_dtriad)))
      shell_dtriad_bak(:) = fstrsolid%shell_dtriad(:)
    endif
    if( associated(fstrsolid%shell_ddrill) ) then
      allocate(shell_ddrill_bak(size(fstrsolid%shell_ddrill)))
      shell_ddrill_bak(:) = fstrsolid%shell_ddrill(:)
    endif
 
    do i=-ntot,ntot
      alpha = 5.d-3*dble(i)/dble(ntot)
      hecmat%X(:) = alpha*fstrsolid%dunode(:)
      fstrsolid%dunode(:) = dunode_bak(:)
      if( allocated(shell_dtriad_bak) ) fstrsolid%shell_dtriad(:) = shell_dtriad_bak(:)
      if( allocated(shell_ddrill_bak) ) fstrsolid%shell_ddrill(:) = shell_ddrill_bak(:)
      call fstr_apply_solution_increment( hecmesh, fstrsolid, hecmat%ndof, hecmat%X )
      call fstr_calc_residual_vector(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam)
      pot(1) = fstr_get_potential(cstep,hecmesh,hecmat,fstrsolid,1)
      pot(2) = fstr_get_potential(cstep,hecmesh,hecmat,fstrsolid,2)
      pot(3) = fstr_get_potential(cstep,hecmesh,hecmat,fstrsolid,3)
      write(imsg,'(I4,5(",",1pE16.9))') i,alpha,alpha*dulen,pot(1:3)
    enddo
 
    do i=1, size(fstrsolid%dunode)
      fstrsolid%dunode(i) = dunode_bak(i)
    enddo
    if( allocated(shell_dtriad_bak) ) then
      fstrsolid%shell_dtriad(:) = shell_dtriad_bak(:)
      deallocate(shell_dtriad_bak)
    endif
    if( allocated(shell_ddrill_bak) ) then
      fstrsolid%shell_ddrill(:) = shell_ddrill_bak(:)
      deallocate(shell_ddrill_bak)
    endif
    hecmat%X(:) = 0.d0
    call fstr_calc_residual_vector_with_x(hecmesh, hecmat, fstrsolid, ctime, tincr, iter, cstep, dtime, fstrparam)
  end subroutine
 
end module m_fstr_residual