mechgauss.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 mmechgauss
  use hecmw_util
  use mmaterial
  use elementinfo, only: fe_mitc3_shell, fe_mitc4_shell, fe_mitc9_shell, numofshellthicknessquadpoints, &
    getshellthicknessquadpoint, getshellthicknessweight
  implicit none
 
  ! ----------------------------------------------------------------------------
  type tgaussstatus
    type(tmaterial), pointer  :: pmaterial => null()    
    real(kind=kreal)          :: strain(6)              
    real(kind=kreal)          :: stress(6)              
    integer, pointer          :: istatus(:) =>null()    
    real(kind=kreal), pointer :: fstatus(:) => null()   
    real(kind=kreal)          :: plstrain               
    real(kind=kreal)          :: strain_bak(6)          
    real(kind=kreal)          :: stress_bak(6)          
    real(kind=kreal)          :: nqm(12)                
    real(kind=kreal)          :: strain_out(6)          
    real(kind=kreal)          :: stress_out(6)          
    real(kind=kreal)          :: strain_energy          
    real(kind=kreal)          :: strain_energy_bak      
    real(kind=kreal)          :: plpotential            
  end type
 
  ! ----------------------------------------------------------------------------
  type telement
    integer                     :: etype
    integer                     :: iset
    real(kind=kreal), pointer   :: equiforces(:) => null()  
    type(tgaussstatus), pointer :: gausses(:) => null()  
    type(tgaussstatus), pointer :: shell_layer_gausses(:) => null()  
    integer(kind=kint)          :: shell_nlayer = 0
    integer(kind=kint)          :: shell_nthick = 0
    real(kind=kreal), pointer   :: aux(:,:) => null()    
    integer                     :: elemact_flag
    real(kind=kreal)            :: elemact_coeff           
    real(kind=kreal)            :: p(1)                  
  end type
 
contains
 
  integer(kind=kint) function fstr_shell_output_point_count(element) result(npoint)
    type(telement), intent(in) :: element
 
    npoint = 0
    if( associated(element%shell_layer_gausses) ) then
      npoint = size(element%shell_layer_gausses)
      return
    endif
    select case(element%etype)
    case(fe_mitc3_shell, fe_mitc4_shell, fe_mitc9_shell)
      if( .not.associated(element%gausses) ) return
      if( element%gausses(1)%pMaterial%nlgeom_flag /= infinitesimal ) return
      if( .not.iselastic(element%gausses(1)%pMaterial%mtype) ) return
      npoint = size(element%gausses)*element%gausses(1)%pMaterial%totallyr &
        *numofshellthicknessquadpoints(element%etype)
    end select
  end function fstr_shell_output_point_count
 
  logical function fstr_uses_plane_stress_j2_shell( etype, nn, material )
    integer(kind=kint), intent(in) :: etype
    integer(kind=kint), intent(in) :: nn
    type(tmaterial), intent(in)    :: material
 
    fstr_uses_plane_stress_j2_shell = .false.
    if( etype /= fe_mitc4_shell .or. nn /= 4 ) return
    if( material%nlgeom_flag /= infinitesimal ) return
    if( .not. iselastoplastic(material%mtype) ) return
    if( getelastictype(material%mtype) /= 0 ) return
    if( getyieldfunction(material%mtype) /= 0 ) return
    if( gethardentype(material%mtype) /= 0 ) return
    if( abs(material%variables(m_plconst2)) > tiny(1.0d0) ) return
    fstr_uses_plane_stress_j2_shell = .true.
  end function fstr_uses_plane_stress_j2_shell
 
  subroutine fstr_init_gauss( gauss )
    use muyield, only: uelastoplasticnumstatus
    type( tgaussstatus ), intent(inout) :: gauss
    integer :: n
    gauss%strain=0.d0; gauss%stress=0.d0
    gauss%strain_bak=0.d0; gauss%stress_bak=0.d0
    gauss%strain_out=0.d0; gauss%stress_out=0.d0
    gauss%plstrain =0.d0
    gauss%plpotential =0.d0
    gauss%nqm =0.d0
    gauss%strain_energy =0.d0
    gauss%strain_energy_bak =0.d0
    if( gauss%pMaterial%mtype==usermaterial ) then
      if( gauss%pMaterial%nfstatus> 0 ) then
        allocate( gauss%fstatus(gauss%pMaterial%nfstatus) )
        gauss%fstatus(:) = 0.d0
      endif
    else if( iselastoplastic(gauss%pMaterial%mtype) ) then
      allocate( gauss%istatus(1) )    ! 0:elastic 1:plastic
      if( getyieldfunction( gauss%pMaterial%mtype )==3 ) then  ! user defined
        n = uelastoplasticnumstatus( gauss%pMaterial%variables )
        if( n>0 ) allocate( gauss%fstatus(n) )
      elseif( iskinematicharden( gauss%pMaterial%mtype ) ) then
        ! (1) trial eq. plastic strain, (2:7) trial back stress, (8:13) converged back stress
        allocate( gauss%fstatus(7+6) )
      else
        allocate( gauss%fstatus(2) )    ! (1) trial eq. plastic strain, (2) unused
      endif
      gauss%istatus = 0
      gauss%fstatus = 0.d0
    else if( isviscoelastic(gauss%pMaterial%mtype) ) then
      n = fetch_tablerow( mc_viscoelastic, gauss%pMaterial%dict )
      if( n>0 ) then
        ! per Prony term: (1:6) converged and (7:12) trial visco stress; last 6: converged deviatoric strain
        allocate( gauss%fstatus(12*n+6) )
        gauss%fstatus = 0.d0
      else
        stop "Viscoelastic properties not defined"
      endif
    else if( gauss%pMaterial%mtype==norton ) then
      allocate( gauss%fstatus(2) )        ! (1) updated eq. stress, (2) converged eq. creep strain
      gauss%fstatus = 0.d0
      gauss%plstrain = 0.d0
    endif
  end subroutine fstr_init_gauss
 
  subroutine fstr_finalize_gauss( gauss )
    type( tgaussstatus ), intent(inout) :: gauss
    if( associated( gauss%istatus ) ) deallocate( gauss%istatus )
    if( associated( gauss%fstatus ) ) deallocate( gauss%fstatus )
  end subroutine
 
  integer(kind=kint) function fstr_shell_num_thickness_points( etype )
    integer(kind=kint), intent(in) :: etype
 
    fstr_shell_num_thickness_points = max(0, numofshellthicknessquadpoints(etype))
  end function fstr_shell_num_thickness_points
 
  subroutine fstr_init_shell_layer_gausses( element, ng, nlayer, nthick )
    type( telement ), intent(inout) :: element
    integer(kind=kint), intent(in)  :: ng, nlayer, nthick
    integer(kind=kint) :: i, nstatus
 
    if( ng <= 0 .or. nlayer <= 0 .or. nthick <= 0 ) return
    if( .not. associated( element%gausses ) ) return
    if( associated( element%shell_layer_gausses ) ) return
 
    nstatus = ng*nlayer*nthick
    allocate( element%shell_layer_gausses( nstatus ) )
    element%shell_nlayer = nlayer
    element%shell_nthick = nthick
 
    do i = 1, nstatus
      element%shell_layer_gausses(i)%pMaterial => element%gausses(1)%pMaterial
      call fstr_init_gauss( element%shell_layer_gausses(i) )
    enddo
  end subroutine fstr_init_shell_layer_gausses
 
  integer(kind=kint) function fstr_shell_layer_gauss_index( element, ig, ilayer, ithick )
    type( telement ), intent(in)    :: element
    integer(kind=kint), intent(in)  :: ig, ilayer, ithick
    integer(kind=kint) :: ngauss
 
    fstr_shell_layer_gauss_index = 0
    if( .not. associated( element%gausses ) ) return
    if( .not. associated( element%shell_layer_gausses ) ) return
    if( element%shell_nlayer <= 0 .or. element%shell_nthick <= 0 ) return
    if( ilayer < 1 .or. ilayer > element%shell_nlayer ) return
    if( ithick < 1 .or. ithick > element%shell_nthick ) return
 
    ngauss = size( element%gausses )
    if( ig < 1 .or. ig > ngauss ) return
    if( size( element%shell_layer_gausses ) < ngauss*element%shell_nlayer*element%shell_nthick ) return
 
    fstr_shell_layer_gauss_index = ((ig-1)*element%shell_nlayer + ilayer-1) &
      *element%shell_nthick + ithick
  end function fstr_shell_layer_gauss_index
 
  real(kind=kreal) function fstr_element_average_plstrain( element )
    type( telement ), intent(in) :: element
    integer(kind=kint) :: i, ig, ilayer, ithick, ishell
    real(kind=kreal) :: weight, total_weight
 
    fstr_element_average_plstrain = 0.0d0
    if( associated(element%shell_layer_gausses) ) then
      if( .not. associated(element%gausses) ) return
      if( element%shell_nlayer <= 0 .or. element%shell_nthick <= 0 ) return
      if( size(element%gausses) <= 0 ) return
      if( .not. associated(element%gausses(1)%pMaterial) ) return
      if( .not. associated(element%gausses(1)%pMaterial%shell_var) ) return
      if( element%shell_nlayer > size(element%gausses(1)%pMaterial%shell_var) ) return
      total_weight = 0.0d0
      do ig = 1, size(element%gausses)
        do ilayer = 1, element%shell_nlayer
          weight = element%gausses(1)%pMaterial%shell_var(ilayer)%weight
          do ithick = 1, element%shell_nthick
            ishell = fstr_shell_layer_gauss_index(element, ig, ilayer, ithick)
            if( ishell <= 0 ) cycle
            fstr_element_average_plstrain = fstr_element_average_plstrain &
              + element%shell_layer_gausses(ishell)%plstrain * weight
            total_weight = total_weight + weight
          enddo
        enddo
      enddo
      if( total_weight > 0.0d0 ) fstr_element_average_plstrain = fstr_element_average_plstrain/total_weight
    else if( associated(element%gausses) ) then
      if( size(element%gausses) <= 0 ) return
      do i = 1, size(element%gausses)
        fstr_element_average_plstrain = fstr_element_average_plstrain+element%gausses(i)%plstrain
      enddo
      fstr_element_average_plstrain = fstr_element_average_plstrain/size(element%gausses)
    endif
  end function fstr_element_average_plstrain
 
  subroutine fstr_shell_thickness_quadrature( etype, ithick, zeta, weight, ierr )
    integer(kind=kint), intent(in)  :: etype, ithick
    real(kind=kreal), intent(out)   :: zeta, weight
    integer(kind=kint), intent(out) :: ierr
 
    ierr = 0
    zeta = 0.0d0
    weight = 0.0d0
 
    if( ithick < 1 .or. ithick > numofshellthicknessquadpoints(etype) ) then
      ierr = 1
      return
    endif
 
    call getshellthicknessquadpoint(etype, ithick, zeta)
    weight = getshellthicknessweight(etype, ithick)
  end subroutine fstr_shell_thickness_quadrature
 
  subroutine fstr_shell_layer_quadrature( element, ilayer, ithick, zeta_layer, weight, ierr )
    type( telement ), intent(in)    :: element
    integer(kind=kint), intent(in)  :: ilayer, ithick
    real(kind=kreal), intent(out)   :: zeta_layer, weight
    integer(kind=kint), intent(out) :: ierr
 
    ierr = 0
    zeta_layer = 0.0d0
    weight = 0.0d0
 
    if( .not. associated( element%gausses ) ) then
      ierr = 1
      return
    endif
 
    call fstr_shell_layer_quadrature_gauss( element%etype, element%gausses(1), ilayer, ithick, &
      zeta_layer, weight, ierr )
  end subroutine fstr_shell_layer_quadrature
 
  subroutine fstr_shell_layer_quadrature_gauss( etype, gauss, ilayer, ithick, zeta_layer, weight, ierr )
    integer(kind=kint), intent(in)  :: etype, ilayer, ithick
    type( tgaussstatus ), intent(in) :: gauss
    real(kind=kreal), intent(out)   :: zeta_layer, weight
    integer(kind=kint), intent(out) :: ierr
    real(kind=kreal) :: zeta
 
    ierr = 0
    zeta_layer = 0.0d0
    weight = 0.0d0
 
    call fstr_shell_thickness_quadrature( etype, ithick, zeta, weight, ierr )
    if( ierr /= 0 ) return
    call fstr_shell_layer_zeta( gauss, ilayer, zeta, zeta_layer, ierr )
  end subroutine fstr_shell_layer_quadrature_gauss
 
  subroutine fstr_shell_layer_zeta( gauss, ilayer, zeta, zeta_layer, ierr )
    type( tgaussstatus ), intent(in) :: gauss
    integer(kind=kint), intent(in)  :: ilayer
    real(kind=kreal), intent(in)    :: zeta
    real(kind=kreal), intent(out)   :: zeta_layer
    integer(kind=kint), intent(out) :: ierr
    integer(kind=kint) :: i
    real(kind=kreal) :: sumlyr
 
    ierr = 0
    zeta_layer = 0.0d0
 
    if( .not. associated( gauss%pMaterial ) ) then
      ierr = 1
      return
    endif
    if( ilayer < 1 .or. ilayer > gauss%pMaterial%totallyr ) then
      ierr = 1
      return
    endif
 
    sumlyr = 0.0d0
    do i = 1, ilayer
      sumlyr = sumlyr + 2.0d0*gauss%pMaterial%shell_var(i)%weight
    enddo
    zeta_layer = -1.0d0 + sumlyr - gauss%pMaterial%shell_var(ilayer)%weight*(1.0d0-zeta)
  end subroutine fstr_shell_layer_zeta
 
  subroutine fstr_finalize_shell_layer_gausses( element )
    type( telement ), intent(inout) :: element
    integer(kind=kint) :: i
 
    if( associated( element%shell_layer_gausses ) ) then
      do i = 1, size( element%shell_layer_gausses )
        call fstr_finalize_gauss( element%shell_layer_gausses(i) )
      enddo
      deallocate( element%shell_layer_gausses )
    endif
    element%shell_nlayer = 0
    element%shell_nthick = 0
  end subroutine fstr_finalize_shell_layer_gausses
 
  subroutine fstr_copy_gauss( gauss1, gauss2 )
    type( tgaussstatus ), intent(in)    :: gauss1
    type( tgaussstatus ), intent(inout) :: gauss2
 
    gauss2%strain     = gauss1%strain
    gauss2%stress     = gauss1%stress
    gauss2%strain_bak = gauss1%strain_bak
    gauss2%stress_bak = gauss1%stress_bak
    gauss2%nqm        = gauss1%nqm
    gauss2%strain_out = gauss1%strain_out
    gauss2%stress_out = gauss1%stress_out
    gauss2%plstrain   = gauss1%plstrain
    gauss2%strain_energy = gauss1%strain_energy
    gauss2%strain_energy_bak = gauss1%strain_energy_bak
    gauss2%plpotential = gauss1%plpotential
 
    if( associated(gauss1%istatus) .and. associated(gauss2%istatus) ) then
      gauss2%istatus   = gauss1%istatus
    end if
    if( associated(gauss1%fstatus) .and. associated(gauss2%fstatus) ) then
      gauss2%fstatus   = gauss1%fstatus
    end if
  end subroutine fstr_copy_gauss
 
  subroutine fstr_copy_shell_layer_gausses( element1, element2 )
    type( telement ), intent(in)    :: element1
    type( telement ), intent(inout) :: element2
    integer(kind=kint) :: i
 
    if( .not. associated( element1%shell_layer_gausses ) ) return
    if( .not. associated( element2%shell_layer_gausses ) ) return
    if( size( element1%shell_layer_gausses ) /= size( element2%shell_layer_gausses ) ) return
 
    do i = 1, size( element1%shell_layer_gausses )
      call fstr_copy_gauss( element1%shell_layer_gausses(i), element2%shell_layer_gausses(i) )
    enddo
  end subroutine fstr_copy_shell_layer_gausses
 
 
end module