fstr_dynamic_nlexplicit.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_dynamic_nlexplicit
  use m_fstr
  use m_static_lib
  use m_dynamic_output
  use m_fstr_eig_setmass
  use m_dynamic_mat_ass_bc_ac
  use m_dynamic_mat_ass_bc
  use m_dynamic_mat_ass_bc_vl
  use m_dynamic_mat_ass_load
  use m_fstr_update
  use m_fstr_restart
  use m_dynamic_mat_ass_couple
  use m_fstr_rcap_io
  use mcontact
 
contains
 
  !C================================================================C
  !C-- subroutine  fstr_solve_LINEAR_DYNAMIC
  !C================================================================C
  subroutine fstr_solve_dynamic_nlexplicit(hecMESH,hecMAT,fstrSOLID,fstrEIG   &
      ,fstrDYN,fstrRESULT,fstrPARAM,infoCTChange &
      ,fstrCPL, restrt_step_num )
    implicit none
    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)                   :: fstrDYN
    type(hecmwst_matrix_lagrange)        :: hecLagMAT
    type(fstr_info_contactchange)        :: infoCTChange
    type(fstr_couple)                    :: fstrCPL !for COUPLE
    type(hecmwst_matrix), pointer :: hecMATmpc
    integer(kind=kint), allocatable :: mark(:)
    integer(kind=kint) :: nnod, ndof, nn, numnp
    integer(kind=kint) :: i, j, ids, ide, kk
    integer(kind=kint) :: kkk0, kkk1
    integer(kind=kint) :: ierror
    integer(kind=kint) :: iiii5, iexit
    integer(kind=kint) :: revocap_flag
    real(kind=kreal), allocatable :: prevb(:)
    real(kind=kreal) :: a1, a2, a3, b1, b2, b3, c1, c2
    real(kind=kreal) :: bsize, res
    real(kind=kreal) :: time_1, time_2
    integer(kind=kint) :: restrt_step_num
    real(kind=kreal), parameter :: pi = 3.14159265358979323846d0
 
    a1 = 0.0d0; a2 = 0.0d0; a3 = 0.0d0; b1 = 0.0d0; b2 = 0.0d0; b3 = 0.0d0
    c1 = 0.0d0; c2 = 0.0d0
 
    call hecmw_mpc_mat_init_explicit(hecmesh, hecmat, hecmatmpc)
 
    hecmat%NDOF=hecmesh%n_dof
    nnod=hecmesh%n_node
    ndof=hecmat%NDOF
    nn=ndof*ndof
 
    if( fstrparam%fg_couple == 1) then
      if( fstrparam%fg_couple_type==5 .or. &
          fstrparam%fg_couple_type==6 ) then
        allocate( prevb(hecmat%NP*ndof)      ,stat=ierror )
        prevb = 0.0d0
        if( ierror /= 0 ) then
          write(idbg,*) 'stop due to allocation error <fstr_solve_NONLINEAR_DYNAMIC, prevB>'
          write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
          call flush(idbg)
          call hecmw_abort( hecmw_comm_get_comm())
        endif
      endif
    endif
 
    fstrsolid%dunode(:) =0.d0
 
    a1 = 1.d0/fstrdyn%t_delta**2
    a2 = 1.d0/(2.d0*fstrdyn%t_delta)
 
    call setmass(fstrsolid,hecmesh,hecmat,fstreig)
    call hecmw_mpc_trans_mass(hecmesh, hecmat, fstreig%mass)
 
    allocate(mark(hecmat%NP * hecmat%NDOF))
    call hecmw_mpc_mark_slave(hecmesh, hecmat, mark)
 
    do j = 1 ,ndof*nnod
      fstrdyn%VEC1(j) = (a1 + a2 *fstrdyn%ray_m) * fstreig%mass(j)
      if(mark(j) == 1) fstrdyn%VEC1(j) = 1.d0
      if(dabs(fstrdyn%VEC1(j)) < 1.0e-20) then
        if( hecmesh%my_rank == 0 ) then
          write(*,*) 'stop due to fstrDYN%VEC(j) = 0 ,  j = ', j
          write(imsg,*) 'stop due to fstrDYN%VEC(j) = 0 ,  j = ', j
        end if
        call hecmw_abort( hecmw_comm_get_comm())
      endif
    end do
 
    deallocate(mark)
 
    !C-- output of initial state
    if( restrt_step_num == 1 ) then
      do j = 1 ,ndof*nnod
        fstrdyn%DISP(j,3) = fstrdyn%DISP(j,1) - fstrdyn%VEL (j,1)/(2.d0*a2)  + fstrdyn%ACC (j,1)/ (2.d0*a1)
        fstrdyn%DISP(j,2) = fstrdyn%DISP(j,1) - fstrdyn%VEL (j,1)/ a2 + fstrdyn%ACC (j,1)/ (2.d0*a1) * 4.d0
      end do
 
      call fstr_dynamic_output(hecmesh, fstrsolid, fstrdyn, fstrparam)
      call dynamic_output_monit(hecmesh, fstrparam, fstrdyn, fstreig, fstrsolid)
    end if
 
    if( associated( fstrsolid%contacts ) )  then
      call initialize_contact_output_vectors(fstrsolid,hecmat)
      call forward_increment_lagrange(1,ndof,fstrdyn%VEC1,hecmesh,fstrsolid,infoctchange,&
        & fstrdyn%DISP(:,2),fstrsolid%ddunode)
    endif
 
    do i= restrt_step_num, fstrdyn%n_step
 
      fstrdyn%i_step = i
      fstrdyn%t_curr = fstrdyn%t_delta * i
 
      !C-- mechanical boundary condition
      call dynamic_mat_ass_load (hecmesh, hecmat, fstrsolid, fstrdyn, fstrparam)
      do j=1, hecmesh%n_node*  hecmesh%n_dof
        hecmat%B(j)=hecmat%B(j)-fstrsolid%QFORCE(j)
      end do
 
      !C ********************************************************************************
      !C for couple analysis
      if( fstrparam%fg_couple == 1 ) then
        if( fstrparam%fg_couple_type==5 .or. &
            fstrparam%fg_couple_type==6 ) then
          do j = 1, hecmat%NP * ndof
            prevb(j) = hecmat%B(j)
          enddo
        endif
      endif
      do
        if( fstrparam%fg_couple == 1 ) then
          if( fstrparam%fg_couple_type==1 .or. &
            fstrparam%fg_couple_type==3 .or. &
            fstrparam%fg_couple_type==5 ) call fstr_rcap_get( fstrcpl )
          if( fstrparam%fg_couple_first /= 0 ) then
            bsize = dfloat( i ) / dfloat( fstrparam%fg_couple_first )
            if( bsize > 1.0 ) bsize = 1.0
            do kkk0 = 1, fstrcpl%coupled_node_n
              kkk1 = 3 * kkk0
              fstrcpl%trac(kkk1-2) = bsize * fstrcpl%trac(kkk1-2)
              fstrcpl%trac(kkk1-1) = bsize * fstrcpl%trac(kkk1-1)
              fstrcpl%trac(kkk1  ) = bsize * fstrcpl%trac(kkk1  )
            enddo
          endif
          if( fstrparam%fg_couple_window > 0 ) then
            j = i - restrt_step_num + 1
            kk = fstrdyn%n_step - restrt_step_num + 1
            bsize = 0.5*(1.0-cos(2.0*pi*dfloat(j)/dfloat(kk)))
            do kkk0 = 1, fstrcpl%coupled_node_n
              kkk1 = 3 * kkk0
              fstrcpl%trac(kkk1-2) = bsize * fstrcpl%trac(kkk1-2)
              fstrcpl%trac(kkk1-1) = bsize * fstrcpl%trac(kkk1-1)
              fstrcpl%trac(kkk1  ) = bsize * fstrcpl%trac(kkk1  )
            enddo
          endif
          call dynamic_mat_ass_couple( hecmesh, hecmat, fstrsolid, fstrcpl )
        endif
        !C ********************************************************************************
 
        call hecmw_mpc_trans_rhs(hecmesh, hecmat, hecmatmpc)
 
        do j = 1 ,ndof*nnod
          hecmatmpc%B(j) = hecmatmpc%B(j) + 2.d0*a1* fstreig%mass(j) * fstrdyn%DISP(j,1)  &
            + (- a1 + a2 * fstrdyn%ray_m) * fstreig%mass(j) * fstrdyn%DISP(j,3)
        end do
 
        !C
        !C-- geometrical boundary condition
 
        call dynamic_explicit_ass_bc(hecmesh, hecmatmpc, fstrsolid, fstrdyn)
        call dynamic_explicit_ass_vl(hecmesh, hecmatmpc, fstrsolid, fstrdyn)
        call dynamic_explicit_ass_ac(hecmesh, hecmatmpc, fstrsolid, fstrdyn)
        !call dynamic_mat_ass_bc   (hecMESH, hecMATmpc, fstrSOLID, fstrDYN, fstrPARAM, hecLagMAT)
        !call dynamic_mat_ass_bc_vl(hecMESH, hecMATmpc, fstrSOLID, fstrDYN, fstrPARAM, hecLagMAT)
        !call dynamic_mat_ass_bc_ac(hecMESH, hecMATmpc, fstrSOLID, fstrDYN, fstrPARAM, hecLagMAT)
 
        ! Finish the calculation
        do j = 1 ,ndof*nnod
          hecmatmpc%X(j) = hecmatmpc%B(j) / fstrdyn%VEC1(j)
          if(dabs(hecmatmpc%X(j)) > 1.0d+5) then
            if( hecmesh%my_rank == 0 ) then
              print *, 'Displacement increment too large, please adjust your step size!',i,hecmatmpc%X(j)
              write(imsg,*) 'Displacement increment too large, please adjust your step size!',i,hecmatmpc%B(j),fstrdyn%VEC1(j)
            end if
            call hecmw_abort( hecmw_comm_get_comm())
          end if
        end do
        call hecmw_mpc_tback_sol(hecmesh, hecmat, hecmatmpc)
 
        !C *****************************************************
        !C for couple analysis
        if( fstrparam%fg_couple == 1 ) then
          if( fstrparam%fg_couple_type>1 ) then
            do j=1, fstrcpl%coupled_node_n
              if( fstrcpl%dof == 3 ) then
                kkk0 = j*3
                kkk1 = fstrcpl%coupled_node(j)*3
 
                fstrcpl%disp (kkk0-2) = hecmat%X(kkk1-2)
                fstrcpl%disp (kkk0-1) = hecmat%X(kkk1-1)
                fstrcpl%disp (kkk0  ) = hecmat%X(kkk1  )
 
                fstrcpl%velo (kkk0-2) = -b1*fstrdyn%ACC(kkk1-2,1) - b2*fstrdyn%VEL(kkk1-2,1) + &
                  b3*( hecmat%X(kkk1-2) - fstrdyn%DISP(kkk1-2,1) )
                fstrcpl%velo (kkk0-1) = -b1*fstrdyn%ACC(kkk1-1,1) - b2*fstrdyn%VEL(kkk1-1,1) + &
                  b3*( hecmat%X(kkk1-1) - fstrdyn%DISP(kkk1-1,1) )
                fstrcpl%velo (kkk0  ) = -b1*fstrdyn%ACC(kkk1,1) - b2*fstrdyn%VEL(kkk1,1) + &
                  b3*( hecmat%X(kkk1) - fstrdyn%DISP(kkk1,1) )
                fstrcpl%accel(kkk0-2) = -a1*fstrdyn%ACC(kkk1-2,1) - a2*fstrdyn%VEL(kkk1-2,1) + &
                  a3*( hecmat%X(kkk1-2) - fstrdyn%DISP(kkk1-2,1) )
                fstrcpl%accel(kkk0-1) = -a1*fstrdyn%ACC(kkk1-1,1) - a2*fstrdyn%VEL(kkk1-1,1) + &
                  a3*( hecmat%X(kkk1-1) - fstrdyn%DISP(kkk1-1,1) )
                fstrcpl%accel(kkk0  ) = -a1*fstrdyn%ACC(kkk1,1) - a2*fstrdyn%VEL(kkk1,1) + &
                  a3*( hecmat%X(kkk1) - fstrdyn%DISP(kkk1,1) )
              else
                kkk0 = j*2
                kkk1 = fstrcpl%coupled_node(j)*2
 
                fstrcpl%disp (kkk0-1) = hecmat%X(kkk1-1)
                fstrcpl%disp (kkk0  ) = hecmat%X(kkk1  )
 
                fstrcpl%velo (kkk0-1) = -b1*fstrdyn%ACC(kkk1-1,1) - b2*fstrdyn%VEL(kkk1-1,1) + &
                  b3*( hecmat%X(kkk1-1) - fstrdyn%DISP(kkk1-1,1) )
                fstrcpl%velo (kkk0  ) = -b1*fstrdyn%ACC(kkk1,1) - b2*fstrdyn%VEL(kkk1,1) + &
                  b3*( hecmat%X(kkk1) - fstrdyn%DISP(kkk1,1) )
                fstrcpl%accel(kkk0-1) = -a1*fstrdyn%ACC(kkk1-1,1) - a2*fstrdyn%VEL(kkk1-1,1) + &
                  a3*( hecmat%X(kkk1-1) - fstrdyn%DISP(kkk1-1,1) )
                fstrcpl%accel(kkk0  ) = -a1*fstrdyn%ACC(kkk1,1) - a2*fstrdyn%VEL(kkk1,1) + &
                  a3*( hecmat%X(kkk1) - fstrdyn%DISP(kkk1,1) )
              endif
            end do
            call fstr_rcap_send( fstrcpl )
          endif
 
          select case ( fstrparam%fg_couple_type )
            case (4)
              call fstr_rcap_get( fstrcpl )
            case (5)
              call fstr_get_convergence( revocap_flag )
              if( revocap_flag==0 ) then
                do j = 1, hecmat%NP * ndof
                  hecmat%B(j) = prevb(j)
                enddo
                cycle
              endif
            case (6)
              call fstr_get_convergence( revocap_flag )
              if( revocap_flag==0 ) then
                do j = 1, hecmat%NP * ndof
                  hecmat%B(j) = prevb(j)
                enddo
                call fstr_rcap_get( fstrcpl )
                cycle
              else
                if( i /= fstrdyn%n_step ) call fstr_rcap_get( fstrcpl )
              endif
          end select
        endif
        exit
      enddo
 
      !C *****************************************************
      !C-- contact corrector
      !C
      do j = 1 ,ndof*nnod
        fstrsolid%unode(j)  = fstrdyn%DISP(j,1)
        fstrsolid%dunode(j)  = hecmat%X(j)-fstrdyn%DISP(j,1)
      enddo
      if( associated( fstrsolid%contacts ) )  then
        !call fstr_scan_contact_state( 1, fstrDYN%t_delta, kcaSLAGRANGE, hecMESH, fstrSOLID, infoCTChange )
        call forward_increment_lagrange(1,ndof,fstrdyn%VEC1,hecmesh,fstrsolid,infoctchange,&
          & fstrdyn%DISP(:,2),fstrsolid%ddunode)
        do j = 1 ,ndof*nnod
          hecmat%X(j)  = hecmat%X(j) + fstrsolid%ddunode(j)
        enddo
      endif
 
      !C-- new displacement, velocity and acceleration
      do j = 1 ,ndof*nnod
        fstrdyn%ACC (j,1) = a1*(hecmat%X(j) - 2.d0*fstrdyn%DISP(j,1) + fstrdyn%DISP(j,3))
        fstrdyn%VEL (j,1) = a2*(hecmat%X(j) - fstrdyn%DISP(j,3))
        fstrsolid%unode(j)  = fstrdyn%DISP(j,1)
        fstrsolid%dunode(j)  = hecmat%X(j)-fstrdyn%DISP(j,1)
        fstrdyn%DISP(j,3) = fstrdyn%DISP(j,1)
        fstrdyn%DISP(j,1) = hecmat%X(j)
        hecmat%X(j)  = fstrsolid%dunode(j)
        fstrdyn%kineticEnergy = fstrdyn%kineticEnergy + 0.5d0*fstreig%mass(j)*fstrdyn%VEL(j,1)*fstrdyn%VEL(j,1)
      end do
 
      ! ----- update strain, stress, and internal force
      call fstr_updatenewton( hecmesh, hecmat, fstrsolid, fstrdyn%t_curr, fstrdyn%t_delta, 0, fstrdyn%strainEnergy )
 
      do j = 1 ,ndof*nnod
        fstrsolid%unode(j) = fstrsolid%unode(j) + fstrsolid%dunode(j)
      end do
      call fstr_updatestate( hecmesh, fstrsolid, fstrdyn%t_delta )
 
      if( fstrdyn%restart_nout > 0 ) then
        if ( mod(i,fstrdyn%restart_nout).eq.0 .or. i.eq.fstrdyn%n_step ) then
          if( associated( fstrsolid%contacts ) )  then
            call fstr_write_restart_dyna_nl(i,hecmesh,fstrsolid,fstrdyn,fstrparam,infoctchange%contactNode_current)
          else
            call fstr_write_restart_dyna_nl(i,hecmesh,fstrsolid,fstrdyn,fstrparam)
          end if
        end if
      end if
      !
      !C-- output new displacement, velocity and acceleration
      call fstr_dynamic_output(hecmesh, fstrsolid, fstrdyn, fstrparam)
      call dynamic_output_monit(hecmesh, fstrparam, fstrdyn, fstreig, fstrsolid)
 
    enddo
 
    if( fstrparam%fg_couple == 1) then
      if( fstrparam%fg_couple_type==5 .or. &
          fstrparam%fg_couple_type==6 ) then
        deallocate( prevb      ,stat=ierror )
        if( ierror /= 0 ) then
          write(idbg,*) 'stop due to deallocation error <fstr_solve_NONLINEAR_DYNAMIC, prevB>'
          write(idbg,*) '  rank = ', hecmesh%my_rank,'  ierror = ',ierror
          call flush(idbg)
          call hecmw_abort( hecmw_comm_get_comm())
        endif
      endif
    endif
 
    call hecmw_mpc_mat_finalize_explicit(hecmesh, hecmat, hecmatmpc)
 
  end subroutine fstr_solve_dynamic_nlexplicit
 
  
  subroutine forward_increment_lagrange(cstep,ndof,mmat,hecMESH,fstrSOLID,infoCTChange,wkarray,uc)
    integer, intent(in)                    :: cstep
    integer, intent(in)                    :: ndof
    real(kind=kreal), intent(in)           :: mmat(:)
    type( hecmwst_local_mesh ), intent(in) :: hecmesh
    type(fstr_solid), intent(inout)        :: fstrSOLID
    type(fstr_info_contactchange)          :: infoCTChange
    real(kind=kreal), intent(out)          :: wkarray(:)
    real(kind=kreal), intent(out)          :: uc(:)
    integer :: i, j, k, m, grpid, slave, nn, iSS, sid, etype, iter
    real(kind=kreal) :: fdum, conv, dlambda, shapefunc(l_max_surface_node), lambda(3)
 
    call fstr_scan_contact_state_exp( cstep, hecmesh, fstrsolid, infoctchange )
    if( .not. infoctchange%active ) return
 
    uc = 0.0d0
 
    iter = 0
    do
      wkarray = 0.0d0
      do i=1,fstrsolid%n_contacts
        do j= 1, size(fstrsolid%contacts(i)%slave)
          if( fstrsolid%contacts(i)%states(j)%state == contactfree ) cycle
          if( fstrsolid%contacts(i)%states(j)%distance>epsilon(1.d0) ) then
            fstrsolid%contacts(i)%states(j)%state = contactfree
            cycle
          endif
          if( iter==0 ) then
            fstrsolid%contacts(i)%states(j)%multiplier(:) =0.d0
            fstrsolid%contacts(i)%states(j)%wkdist =0.d0
            cycle
          endif
          slave = fstrsolid%contacts(i)%slave(j)
 
          sid = fstrsolid%contacts(i)%states(j)%surface
          nn = size( fstrsolid%contacts(i)%master(sid)%nodes )
          etype = fstrsolid%contacts(i)%master(sid)%etype
          call getshapefunc( etype, fstrsolid%contacts(i)%states(j)%lpos(:), shapefunc )
          wkarray( slave ) = -fstrsolid%contacts(i)%states(j)%multiplier(1)
          do k=1,nn
            iss = fstrsolid%contacts(i)%master(sid)%nodes(k)
            wkarray( iss ) = wkarray( iss ) + shapefunc(k) * fstrsolid%contacts(i)%states(j)%multiplier(1)
          enddo
        enddo
      enddo
 
      if(iter > 0)then
        do i=1,fstrsolid%n_contacts
          do j= 1, size(fstrsolid%contacts(i)%slave)
            if( fstrsolid%contacts(i)%states(j)%state == contactfree ) cycle
            slave = fstrsolid%contacts(i)%slave(j)
            sid = fstrsolid%contacts(i)%states(j)%surface
            nn = size( fstrsolid%contacts(i)%master(sid)%nodes )
            etype = fstrsolid%contacts(i)%master(sid)%etype
            call getshapefunc( etype, fstrsolid%contacts(i)%states(j)%lpos(:), shapefunc )
            fstrsolid%contacts(i)%states(j)%wkdist = -wkarray( slave )/mmat( (slave-1)*ndof+1 )
            do k=1,nn
              iss = fstrsolid%contacts(i)%master(sid)%nodes(k)
              fstrsolid%contacts(i)%states(j)%wkdist = fstrsolid%contacts(i)%states(j)%wkdist  &
                   + shapefunc(k) * wkarray(iss) / mmat( (iss-1)*ndof+1 )
            enddo
          enddo
        enddo
      endif
 
      conv = 0.d0
      wkarray = 0.d0
      do i=1,fstrsolid%n_contacts
        do j= 1, size(fstrsolid%contacts(i)%slave)
          if( fstrsolid%contacts(i)%states(j)%state == contactfree ) cycle
          slave = fstrsolid%contacts(i)%slave(j)
          sid = fstrsolid%contacts(i)%states(j)%surface
          nn = size( fstrsolid%contacts(i)%master(sid)%nodes )
          etype = fstrsolid%contacts(i)%master(sid)%etype
          call getshapefunc( etype, fstrsolid%contacts(i)%states(j)%lpos(:), shapefunc )
          fdum = 1.d0/mmat( (slave-1)*ndof+1 )
          do k=1,nn
            iss = fstrsolid%contacts(i)%master(sid)%nodes(k)
            fdum = fdum + shapefunc(k)*shapefunc(k)/mmat( (iss-1)*ndof+1 )
          enddo
          dlambda= (fstrsolid%contacts(i)%states(j)%distance-fstrsolid%contacts(i)%states(j)%wkdist) /fdum
          conv = conv + dlambda*dlambda;
          fstrsolid%contacts(i)%states(j)%multiplier(1) = fstrsolid%contacts(i)%states(j)%multiplier(1) + dlambda
          if( fstrsolid%contacts(i)%fcoeff>0.d0 ) then
            if( fstrsolid%contacts(i)%states(j)%state == contactslip ) then
              fstrsolid%contacts(i)%states(j)%multiplier(2) =             &
              fstrsolid%contacts(i)%fcoeff * fstrsolid%contacts(i)%states(j)%multiplier(1)
            else    ! stick
              !      fstrSOLID%contacts(i)%states(j)%multiplier(2) =
            endif
          endif
          lambda = fstrsolid%contacts(i)%states(j)%multiplier(1)* fstrsolid%contacts(i)%states(j)%direction
          wkarray((slave-1)*ndof+1:(slave-1)*ndof+3) = lambda(:)
          do k=1,nn
            iss = fstrsolid%contacts(i)%master(sid)%nodes(k)
            wkarray((iss-1)*ndof+1:(iss-1)*ndof+3) = wkarray((iss-1)*ndof+1:(iss-1)*ndof+3) -lambda(:)*shapefunc(k)
          enddo
        enddo
      enddo
      if( dsqrt(conv)<1.d-8 ) exit
      iter = iter+1
    enddo
 
    do i=1,hecmesh%n_node*ndof
      uc(i) = wkarray(i)/mmat(i)
    enddo
  end subroutine forward_increment_lagrange
 
end module fstr_dynamic_nlexplicit