fstr_Update.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_update
  use m_fstr
  implicit none
 
  private :: update_abort
 
contains
 
  !=====================================================================*
  !  UPDATE_C3
  subroutine fstr_updatenewton ( hecMESH, hecMAT, fstrSOLID, time, tincr,iter, strainEnergy)
    !=====================================================================*
    use m_static_lib
    use m_elemact
    use m_fstr_nodalkinematics, only: fstr_ensure_finite_rotation_state, fstr_get_shell_trial_directors, &
      fstr_get_shell_reference_directors
    use m_fstr_finiterotationkinematics, only: fstr_uses_finite_rotation_kinematics
 
    type (hecmwST_matrix)       :: hecMAT
    type (hecmwST_local_mesh)   :: hecMESH
    type (fstr_solid)           :: fstrSOLID
    real(kind=kreal),intent(in) :: time      
    real(kind=kreal),intent(in) :: tincr     
    integer, intent(in)         :: iter
 
    integer(kind=kint) :: nodLOCAL(fstrSOLID%max_ncon)
    real(kind=kreal)   :: ecoord(3, fstrsolid%max_ncon)
    real(kind=kreal)   :: thick, thick0(6)
    integer(kind=kint) :: ndof, itype, is, iE, ic_type, nn, icel, iiS, i, j
 
    real(kind=kreal)   :: total_disp(6, fstrsolid%max_ncon), du(6, fstrsolid%max_ncon), ddu(6, fstrsolid%max_ncon)
    real(kind=kreal)   :: shell_director(3, fstrsolid%max_ncon)
    real(kind=kreal)   :: shell_ref_director(3, fstrsolid%max_ncon)
    real(kind=kreal)   :: tt(fstrsolid%max_ncon), tt0(fstrsolid%max_ncon), ttn(fstrsolid%max_ncon)
    real(kind=kreal)   :: qf(fstrsolid%max_ncon*6), coords(3, 3)
    integer            :: isect, ihead, cdsys_ID
    integer            :: ndim, initt
 
    real(kind=kreal), optional :: strainenergy
    real(kind=kreal) :: tmp, lambda(1), ddlambda(1)
    real(kind=kreal)   :: ddaux(3,3)
 
    ndof = hecmat%NDOF
    fstrsolid%QFORCE=0.0d0
    call fstr_ensure_finite_rotation_state( hecmesh, fstrsolid, ndof )
 
    tt0 = 0.d0
    ttn = 0.d0
    tt = 0.d0
 
    ! if initial temperature exists
    initt = 0
    if( associated(g_initialcnd) ) then
        do j=1,size(g_initialcnd)
          if( g_initialcnd(j)%cond_name=="temperature" ) then
            initt=j
            exit
          endif
        end do
    endif
 
    ! --------------------------------------------------------------------
    !      updated
    !   1. stress and strain  : ep^(k) = ep^(k-1)+dep^(k)
    !                           sgm^(k) = sgm^(k-1)+dsgm^(k)
    !   2. Internal Force     : Q^(k-1) ( u^(k-1) )
    ! --------------------------------------------------------------------
    ! ----------------------------------------------------------------------------------
    !      calculate the Strain and Stress and Internal Force ( Equivalent Nodal Force )
    ! ----------------------------------------------------------------------------------
 
    do itype = 1, hecmesh%n_elem_type
      is = hecmesh%elem_type_index(itype-1)+1
      ie = hecmesh%elem_type_index(itype  )
      ic_type= hecmesh%elem_type_item(itype)
      if (hecmw_is_etype_link(ic_type)) cycle
      if (hecmw_is_etype_patch(ic_type)) cycle
 
      !element loop
      !$omp parallel default(none), &
        !$omp&  private(icel,iiS,j,nn,nodLOCAL,i,ecoord,ddu,du,total_disp,shell_director,shell_ref_director, &
        !$omp&  cdsys_ID,coords,thick,qf,isect,ihead,tmp,ndim,ddaux,thick0, &
        !$omp&  lambda,ddlambda), &
        !$omp&  shared(iS,iE,hecMESH,fstrSOLID,ndof,hecMAT,ic_type,fstrPR, &
        !$omp&         strainEnergy,iter,time,tincr,initt,g_InitialCnd), &
        !$omp&  firstprivate(tt0,ttn,tt)
      !$omp do
      do icel = is, ie
 
        ! ----- nodal coordinate, displacement and temperature
        iis = hecmesh%elem_node_index(icel-1)
        nn = hecmesh%elem_node_index(icel)-iis
        !if( nn>150 ) stop "elemental nodes > 150!"
 
        thick = 0.d0
        do j = 1, size(fstrsolid%elements(icel)%gausses)
          thick0(1:6) = fstrsolid%elements(icel)%gausses(j)%stress(1:6)
          thick = thick + dsqrt(dot_product(thick0(1:6),thick0(1:6)))
        enddo
        if( thick < 1.d-10 ) then
          do j = 1, size(fstrsolid%elements(icel)%gausses)
            if( associated(fstrsolid%elements(icel)%gausses(j)%fstatus) ) &
              &  fstrsolid%elements(icel)%gausses(j)%fstatus = 0.d0
          enddo
        end if
 
        do j = 1, nn
          nodlocal(j) = hecmesh%elem_node_item (iis+j)
          do i = 1, 3
            ecoord(i,j) = hecmesh%node(3*nodlocal(j)+i-3)
          enddo
          do i = 1, ndof
            ddu(i,j) = hecmat%X(ndof*nodlocal(j)+i-ndof)
            du(i,j)  = fstrsolid%dunode(ndof*nodlocal(j)+i-ndof)
            total_disp(i,j) = fstrsolid%unode(ndof*nodlocal(j)+i-ndof)
          enddo
 
          if( fstrsolid%TEMP_ngrp_tot > 0 .or. fstrsolid%TEMP_irres > 0 ) then
            if( iselastoplastic(fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype) .or. &
                fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype == norton ) then
              tt0(j)=fstrsolid%last_temp( nodlocal(j) )
            else
              tt0(j) = 0.d0
              if( hecmesh%hecmw_flag_initcon == 1 ) tt0(j) = hecmesh%node_init_val_item(nodlocal(j))
              if( initt>0 ) tt0(j) = g_initialcnd(initt)%realval(nodlocal(j))
            endif
            ttn(j) = fstrsolid%last_temp( nodlocal(j) )
            tt(j)  = fstrsolid%temperature( nodlocal(j) )
          endif
        enddo
 
        isect = hecmesh%section_ID(icel)
        ihead = hecmesh%section%sect_R_index(isect-1)
        thick = hecmesh%section%sect_R_item(ihead+1)
        cdsys_id = hecmesh%section%sect_orien_ID(isect)
        if( cdsys_id > 0 ) call get_coordsys(cdsys_id, hecmesh, fstrsolid, coords)
 
        ! ===== calculate the Internal Force
        if( ic_type == 241 .or. ic_type == 242 .or. ic_type == 231 .or. ic_type == 232 .or. ic_type == 2322 ) then
          call update_c2( ic_type,nn,ecoord(1:3,1:nn),fstrsolid%elements(icel)%gausses(:), &
            thick,fstrsolid%elements(icel)%iset,                             &
            total_disp(1:2,1:nn), ddu(1:2,1:nn), qf(1:nn*ndof),              &
            tt(1:nn), tt0(1:nn), ttn(1:nn)  )
 
        else if( ic_type == 301 ) then
          call update_c1( ic_type,nn,ecoord(:,1:nn), thick, total_disp(1:3,1:nn), du(1:3,1:nn), &
            qf(1:nn*ndof),fstrsolid%elements(icel)%gausses(:) )
 
        else if( ic_type == 361 ) then
          if( fstrsolid%sections(isect)%elemopt361 == kel361fi ) then ! full integration element
            call update_c3( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
              qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn)  )
          else if( fstrsolid%sections(isect)%elemopt361 == kel361bbar ) then ! B-bar element
            call update_c3d8bbar( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords,    &
              qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn)  )
          else if( fstrsolid%sections(isect)%elemopt361 == kel361ic ) then ! incompatible element
            call update_c3d8ic( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), ddu(1:3,1:nn), cdsys_id, coords,&
              qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, &
              fstrsolid%elements(icel)%aux, ddaux(1:3,1:3), tt(1:nn), tt0(1:nn), ttn(1:nn) )
            fstrsolid%elements(icel)%aux(1:3,1:3) = fstrsolid%elements(icel)%aux(1:3,1:3) + ddaux(1:3,1:3)
          else if( fstrsolid%sections(isect)%elemopt361 == kel361fbar ) then ! F-bar element
            call update_c3d8fbar( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords,    &
              qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn)  )
          else if( fstrsolid%sections(isect)%elemopt361 == kel361up ) then ! UP element
              lambda(1) = -0.5d0*fstrsolid%elements(icel)%p(1)
              CALL update_c3_up                                                                      &
                   ( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), ddu(1:3,1:nn), &
                     cdsys_id, coords, qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:),           &
                     iter, time, tincr, 1, lambda, ddlambda, tt(1:nn), tt0(1:nn) )
              lambda(1) = lambda(1) + ddlambda(1)
              fstrsolid%elements(icel)%p(1) = -2.0d0*lambda(1)
          endif
        else if (ic_type == 341 .or. ic_type == 351 .or. ic_type == 342 .or. ic_type == 352 .or. ic_type == 362 ) then
          if( ic_type==341 .and. fstrsolid%sections(isect)%elemopt341 == kel341sesns ) cycle ! skip smoothed fem
          call update_c3( ic_type, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
            qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), iter, time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn)  )
 
        else if( ic_type == 511) then
          call update_connector( ic_type,nn,ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), &
            qf(1:nn*ndof),fstrsolid%elements(icel)%gausses(:) )
  
        else if( ic_type == 611) then
          if( fstrpr%nlgeom ) call update_abort( ic_type, 2 )
          CALL updatest_beam(ic_type, nn, ecoord, total_disp(1:6,1:nn), du(1:6,1:nn), &
                 &   hecmesh%section%sect_R_item(ihead+1:), fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof))
 
        else if( ic_type == 641 ) then
          if( fstrpr%nlgeom ) call update_abort( ic_type, 2 )
          call updatest_beam_641(ic_type, nn, ecoord, total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
            &    fstrsolid%elements(icel)%gausses(:), hecmesh%section%sect_R_item(ihead+1:), qf(1:nn*ndof))
 
        else if( ( ic_type == 741 ) .or. ( ic_type == 743 ) .or. ( ic_type == 731 ) ) then
          if( fstr_uses_finite_rotation_kinematics( ic_type, nn, &
              fstrsolid%elements(icel)%gausses(1)%pMaterial ) ) then
            call fstr_get_shell_trial_directors( fstrsolid, thick, nn, nodlocal(1:nn), shell_director(1:3,1:nn) )
            call fstr_get_shell_reference_directors( fstrsolid, thick, nn, nodlocal(1:nn), &
              shell_ref_director(1:3,1:nn) )
            call updatest_shell_mitc(ic_type, nn, ndof, ecoord(1:3, 1:nn), total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
              &              fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 0, &
              &              element=fstrsolid%elements(icel), nddirector=shell_director(1:3,1:nn), &
              &              ndrefdirector=shell_ref_director(1:3,1:nn))
          else
            if( fstrpr%nlgeom ) call update_abort( ic_type, 2 )
            call updatest_shell_mitc(ic_type, nn, ndof, ecoord(1:3, 1:nn), total_disp(1:ndof,1:nn), &
              &              du(1:ndof,1:nn), fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 0)
          endif
 
        else if( ic_type == 761 ) then   !for shell-solid mixed analysis
          if( fstrpr%nlgeom ) call update_abort( ic_type, 2 )
          call updatest_shell_mitc33(731, 3, 6, ecoord(1:3, 1:3), total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
            &              fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 2)
 
        else if( ic_type == 781 ) then   !for shell-solid mixed analysis
          if( fstrpr%nlgeom ) call update_abort( ic_type, 2 )
          call updatest_shell_mitc33(741, 4, 6, ecoord(1:3, 1:4), total_disp(1:ndof,1:nn), du(1:ndof,1:nn), &
            &              fstrsolid%elements(icel)%gausses(:), qf(1:nn*ndof), thick, 1)
 
        else if ( ic_type == 3414 ) then
          if(fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype /= incomp_newtonian) &
            &  call update_abort( ic_type, 3, fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype )
          call update_c3_vp                                                       &
            ( ic_type, nn, ecoord(:,1:nn), total_disp(1:4,1:nn), du(1:4,1:nn), &
            fstrsolid%elements(icel)%gausses(:) )
          qf = 0.0d0
 
        else if ( ic_type == 881 .or. ic_type == 891 ) then  !for selective es/ns smoothed fem
          call update_c3_sesns( ic_type, nn, nodlocal, ecoord(:,1:nn), total_disp(1:3,1:nn), du(1:3,1:nn), cdsys_id, coords, &
            qf(1:nn*ndof), fstrsolid%elements(icel)%gausses(:), time, tincr, tt(1:nn), tt0(1:nn), ttn(1:nn)  )
 
        else
          write(*, *) '###ERROR### : Element type not supported for nonlinear static analysis'
          write(*, *) ' ic_type = ', ic_type
          call hecmw_abort(hecmw_comm_get_comm())
 
        endif
 
        ! elemact element
        if( fstrsolid%elements(icel)%elemact_flag == kelact_inactive ) then
          call update_dummy( ndof, nn, ecoord(:,1:nn), total_disp(1:3,1:nn), &
            &  du(1:3,1:nn), qf(1:nn*ndof), fstrsolid%elements(icel) )
          !qf(:) = fstrSOLID%elements(icel)%elemact_coeff*qf(:)
        end if
 
        ! ----- calculate the global internal force ( Q(u_{n+1}^{k-1}) )
        do j = 1, nn
          do i = 1, ndof
            !$omp atomic
            fstrsolid%QFORCE(ndof*(nodlocal(j)-1)+i) = fstrsolid%QFORCE(ndof*(nodlocal(j)-1)+i)+qf(ndof*(j-1)+i)
          enddo
        enddo
 
        ! ----- calculate strain energy
        if(present(strainenergy))then
          ndim = getspacedimension( fstrsolid%elements(icel)%etype )
          do j = 1, nn
            do i = 1, ndim
              tmp = 0.5d0*( fstrsolid%elements(icel)%equiForces(ndim*(j-1)+i)+qf(ndim*(j-1)+i) )*ddu(i,j)
              !$omp atomic
              strainenergy = strainenergy+tmp
              fstrsolid%elements(icel)%equiForces(ndim*(j-1)+i) = qf(ndim*(j-1)+i)
            enddo
          enddo
        endif
 
      enddo ! icel
      !$omp end do
      !$omp end parallel
    enddo   ! itype
 
    !C
    !C Update for fstrSOLID%QFORCE
    !C
    call hecmw_update_r(hecmesh,fstrsolid%QFORCE,hecmesh%n_node, ndof)
  end subroutine fstr_updatenewton
 
 
  subroutine fstr_updatestate( hecMESH, fstrSOLID, tincr)
    use m_fstr
    use m_static_lib
    use m_elastoplastic
    use mcreep
    use mviscoelastic
    type(hecmwst_local_mesh) :: hecmesh
    type(fstr_solid) :: fstrSOLID
    real(kind=kreal) :: tincr
    integer(kind=kint) :: itype, is, iE, ic_type, icel, ngauss, i
    integer(kind=kint) :: ishell
 
    if( associated( fstrsolid%temperature ) ) then
      do i = 1, hecmesh%n_node
        fstrsolid%last_temp(i) = fstrsolid%temperature(i)
      end do
    endif
 
    do itype = 1, hecmesh%n_elem_type
      is = hecmesh%elem_type_index(itype-1) + 1
      ie = hecmesh%elem_type_index(itype  )
      ic_type= hecmesh%elem_type_item(itype)
      if( ic_type == 301 ) ic_type = 111
      if( hecmw_is_etype_link(ic_type) ) cycle
      if( hecmw_is_etype_patch(ic_type) ) cycle
 
      ngauss = numofquadpoints( ic_type )
      do icel = is, ie
        if( iselastoplastic( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype ) ) then
          do i = 1, ngauss
            call updateepstate( fstrsolid%elements(icel)%gausses(i) )
          enddo
        elseif( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype == norton ) then
          if( tincr>0.d0 ) then
            do i = 1, ngauss
              call updateviscostate( fstrsolid%elements(icel)%gausses(i) )
            enddo
          endif
        elseif( isviscoelastic( fstrsolid%elements(icel)%gausses(1)%pMaterial%mtype ) ) then
          if( tincr > 0.d0 ) then
            do i = 1, ngauss
              call updateviscoelasticstate( fstrsolid%elements(icel)%gausses(i) )
            enddo
          endif
        endif
 
        do i = 1, ngauss
          fstrsolid%elements(icel)%gausses(i)%strain_bak = fstrsolid%elements(icel)%gausses(i)%strain
          fstrsolid%elements(icel)%gausses(i)%stress_bak = fstrsolid%elements(icel)%gausses(i)%stress
          fstrsolid%elements(icel)%gausses(i)%strain_energy_bak = fstrsolid%elements(icel)%gausses(i)%strain_energy
        enddo
 
        if( associated( fstrsolid%elements(icel)%shell_layer_gausses ) ) then
          do ishell = 1, size( fstrsolid%elements(icel)%shell_layer_gausses )
            fstrsolid%elements(icel)%shell_layer_gausses(ishell)%strain_bak = &
              fstrsolid%elements(icel)%shell_layer_gausses(ishell)%strain
            fstrsolid%elements(icel)%shell_layer_gausses(ishell)%stress_bak = &
              fstrsolid%elements(icel)%shell_layer_gausses(ishell)%stress
            fstrsolid%elements(icel)%shell_layer_gausses(ishell)%strain_energy_bak = &
              fstrsolid%elements(icel)%shell_layer_gausses(ishell)%strain_energy
          enddo
        endif
      enddo
    enddo
 
    do i = 1, hecmesh%n_node
      fstrsolid%QFORCE_bak(i) = fstrsolid%QFORCE(i)
    end do
 
  end subroutine fstr_updatestate
 
  subroutine update_abort( ic_type, flag, mtype )
    integer(kind=kint), intent(in) :: ic_type
    integer(kind=kint), intent(in) :: flag
    integer(kind=kint), intent(in), optional :: mtype
 
    if( flag == 1 ) then
      write(*,*) '###ERROR### : Element type not supported for static analysis'
    else if( flag == 2 ) then
      write(*,*) '###ERROR### : Element type not supported for nonlinear static analysis'
    else if( flag == 3 ) then
      write(*,*) '###ERROR### : This element is not supported for this material'
    endif
    write(*,*) ' ic_type = ', ic_type
    if( present(mtype) ) write(*,*) ' mtype = ', mtype
    call hecmw_abort(hecmw_comm_get_comm())
  end subroutine
 
end module m_fstr_update