calMatMatrix.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_matmatrix
  use hecmw_util
  use mmaterial
  use mmechgauss
  use m_elasticlinear
  use mhyperelastic
  use m_elastoplastic
  use mviscoelastic
  use mcreep
  use muelastic
  use mumat
 
  implicit none
 
contains
 
  integer function getnlgeomflag( gauss )
    type( tgaussstatus ), intent(in) :: gauss
    getnlgeomflag = gauss%pMaterial%nlgeom_flag
  end function
 
  subroutine matlmatrix( gauss, sectType, matrix, time, dtime, cdsys, temperature, isEp, hdflag )
    type( tgaussstatus ), intent(in) :: gauss
    integer, intent(in)              :: sectType
    real(kind=kreal), intent(out)    :: matrix(:,:)
    real(kind=kreal), intent(in)     :: time
    real(kind=kreal), intent(in)     :: dtime
    real(kind=kreal), intent(in)     :: cdsys(3,3)
    real(kind=kreal), intent(in)     :: temperature
    integer(kind=kint), intent(in), optional :: isEp
    integer(kind=kint), intent(in), optional :: hdflag
 
    integer :: i
    integer :: flag, hdflag_in
    real(kind=kreal)            :: cijkl(3,3,3,3)
    type( tmaterial ), pointer  :: matl
    matl=>gauss%pMaterial
 
    flag = 0
    if( present(isep) )then
      if( isep == 1 )flag = 1
    endif
 
    hdflag_in = 0
    if( present(hdflag) ) hdflag_in = hdflag
 
    if( matl%mtype==userelastic ) then
      call uelasticmatrix( matl%variables(101:), gauss%strain, matrix )
    elseif( isviscoelastic(matl%mtype) ) then
      call calviscoelasticmatrix( matl, secttype, dtime, matrix, temperature )
    elseif( iselastic(matl%mtype) .or. flag==1 ) then
      if(flag==1)then
        i = getelastictype(elastic)
      else
        i = getelastictype(gauss%pMaterial%mtype)
      endif
 
      if( i==0 ) then
        call calelasticmatrix( matl, secttype, matrix, temperature, hdflag=hdflag_in )
      elseif(  i==1 ) then
        call calelasticmatrix_ortho( gauss%pMaterial, secttype, cdsys, matrix, temperature )
      else
        print *, "Elasticity type", matl%mtype, "not supported"
        stop
      endif
 
    elseif( matl%mtype==neohooke .or. matl%mtype==mooneyrivlin ) then
      call calelasticmooneyrivlin( matl, secttype, cijkl, gauss%strain, hdflag=hdflag_in )
      call mat_c2d( cijkl, matrix, secttype )
    elseif( matl%mtype==arrudaboyce )  then
      call calelasticarrudaboyce( matl, secttype, cijkl, gauss%strain )
      call mat_c2d( cijkl, matrix, secttype )
    elseif( matl%mtype==mooneyrivlin_aniso ) then
      call calelasticmooneyrivlinaniso( matl, secttype, cijkl, gauss%strain, cdsys, hdflag=hdflag_in )
      call mat_c2d( cijkl, matrix, secttype )
    elseif( matl%mtype==userhyperelastic )  then
      call uelasticmatrix( matl%variables(101:), gauss%strain, matrix )
    elseif( iselastoplastic(matl%mtype) )  then
      call calelastoplasticmatrix( matl, secttype, gauss%stress,  &
        gauss%istatus(1), gauss%fstatus, gauss%plstrain, matrix, temperature, hdflag=hdflag_in  )
    elseif( matl%mtype==usermaterial ) then
      call umatlmatrix( matl%name, matl%variables(101:), gauss%strain,  &
        gauss%stress, gauss%fstatus, matrix, dtime, time )
    elseif( matl%mtype==norton ) then
      call iso_creep( matl, secttype, gauss%stress, gauss%strain, gauss%fstatus,  &
        gauss%plstrain, dtime, time, matrix, temperature, hdflag=hdflag_in )
    else
      stop "Material type not supported!"
    endif
 
  end subroutine
 
  !
  subroutine stressupdate( gauss, sectType, strain, stress, cdsys, time, dtime, temp, tempn, hdflag  )
    type( tgaussstatus ), intent(inout) :: gauss
    integer, intent(in)                 :: sectType
    real(kind=kreal), intent(in)        :: strain(6)
    real(kind=kreal), intent(out)       :: stress(6)
    real(kind=kreal), intent(in)        :: cdsys(3,3)
    real(kind=kreal), intent(in), optional  :: time
    real(kind=kreal), intent(in), optional  :: dtime
    real(kind=kreal), intent(in)        :: temp
    real(kind=kreal), intent(in)        :: tempn
    integer(kind=kint), intent(in)      :: hdflag
 
    if( gauss%pMaterial%mtype==neohooke .or. gauss%pMaterial%mtype==mooneyrivlin ) then
      call calupdateelasticmooneyrivlin( gauss%pMaterial, secttype, strain, stress, gauss%strain_energy, hdflag=hdflag )
    elseif( gauss%pMaterial%mtype==arrudaboyce ) then ! Arruda-Boyce Hyperelastic material
      call calupdateelasticarrudaboyce( gauss%pMaterial, secttype, strain, stress )
    elseif( gauss%pMaterial%mtype==mooneyrivlin_aniso ) then
      call calupdateelasticmooneyrivlinaniso( gauss%pMaterial, secttype, strain, stress, cdsys, hdflag=hdflag )
    elseif( gauss%pMaterial%mtype==userhyperelastic .or. gauss%pMaterial%mtype==userelastic ) then ! user-defined
      call uelasticupdate( gauss%pMaterial%variables(101:), strain, stress )
    elseif( isviscoelastic( gauss%pMaterial%mtype) ) then
      if( .not. present(dtime) ) stop "error in viscoelastic update!"
      call updateviscoelastic( gauss%pMaterial, secttype, strain, stress, gauss%fstatus, dtime, temp, tempn )
    elseif ( gauss%pMaterial%mtype==norton ) then
      if( .not. present(dtime)  ) stop "error in viscoelastic update!"
      call update_iso_creep( gauss%pMaterial, secttype, strain, stress, gauss%fstatus, &
        &  gauss%plstrain, dtime, time, temp, hdflag=hdflag )
    elseif ( gauss%pMaterial%mtype==usermaterial)  then ! user-defined
      call uupdate(  gauss%pMaterial%name, gauss%pMaterial%variables(101:),   &
        strain, stress, gauss%fstatus, dtime, time )
    end if
 
  end subroutine stressupdate
 
  subroutine mat_c2d( cijkl, dij, itype )
    real(kind=kreal),   intent(in)  :: cijkl(3,3,3,3)
    real(kind=kreal),   intent(out) :: dij(6,6)
    integer,            intent(in)  :: itype
 
    dij(:,:) = 0.d0
    select case( itype )
      case( d3 )
        dij(1,1) = cijkl(1,1,1,1)   ! ---
        dij(1,2) = cijkl(1,1,2,2)
        dij(1,3) = cijkl(1,1,3,3)
        dij(1,4) = cijkl(1,1,1,2)
        dij(1,5) = cijkl(1,1,2,3)
        dij(1,6) = cijkl(1,1,3,1)
        dij(2,1) = cijkl(2,2,1,1)   ! ---
        dij(2,2) = cijkl(2,2,2,2)
        dij(2,3) = cijkl(2,2,3,3)
        dij(2,4) = cijkl(2,2,1,2)
        dij(2,5) = cijkl(2,2,2,3)
        dij(2,6) = cijkl(2,2,3,1)
        dij(3,1) = cijkl(3,3,1,1)   ! ---
        dij(3,2) = cijkl(3,3,2,2)
        dij(3,3) = cijkl(3,3,3,3)
        dij(3,4) = cijkl(3,3,1,2)
        dij(3,5) = cijkl(3,3,2,3)
        dij(3,6) = cijkl(3,3,3,1)
        dij(4,1) = cijkl(1,2,1,1)   ! ---
        dij(4,2) = cijkl(1,2,2,2)
        dij(4,3) = cijkl(1,2,3,3)
        dij(4,4) = cijkl(1,2,1,2)
        dij(4,5) = cijkl(1,2,2,3)
        dij(4,6) = cijkl(1,2,3,1)
        dij(5,1) = cijkl(2,3,1,1)   ! ---
        dij(5,2) = cijkl(2,3,2,2)
        dij(5,3) = cijkl(2,3,3,3)
        dij(5,4) = cijkl(2,3,1,2)
        dij(5,5) = cijkl(2,3,2,3)
        dij(5,6) = cijkl(2,3,3,1)
        dij(6,1) = cijkl(3,1,1,1)   ! ---
        dij(6,2) = cijkl(3,1,2,2)
        dij(6,3) = cijkl(3,1,3,3)
        dij(6,4) = cijkl(3,1,1,2)
        dij(6,5) = cijkl(3,1,2,3)
        dij(6,6) = cijkl(3,1,3,1)
        !
      case( planestress, planestrain )
        dij(1,1) = cijkl(1,1,1,1)   ! ---
        dij(1,2) = cijkl(1,1,2,2)
        dij(1,3) = cijkl(1,1,1,2)
        dij(2,1) = cijkl(2,2,1,1)   ! ---
        dij(2,2) = cijkl(2,2,2,2)
        dij(2,3) = cijkl(2,2,1,2)
        dij(3,1) = cijkl(1,2,1,1)   ! ---
        dij(3,2) = cijkl(1,2,2,2)
        dij(3,3) = cijkl(1,2,1,2)
      case( axissymetric )
        dij(1,1) = cijkl(1,1,1,1)
        dij(1,2) = cijkl(1,1,2,2)
        dij(1,3) = cijkl(1,1,1,2)
        dij(1,4) = cijkl(1,1,3,3)
        dij(2,1) = cijkl(2,2,1,1)
        dij(2,2) = cijkl(2,2,2,2)
        dij(2,3) = cijkl(2,2,1,2)
        dij(2,4) = cijkl(2,2,3,3)
        dij(3,1) = cijkl(1,2,1,1)
        dij(3,2) = cijkl(1,2,2,2)
        dij(3,3) = cijkl(1,2,1,2)
        dij(3,4) = cijkl(1,2,3,3)
        dij(4,1) = cijkl(3,3,1,1)
        dij(4,2) = cijkl(3,3,2,2)
        dij(4,3) = cijkl(3,3,1,2)
        dij(4,4) = cijkl(3,3,3,3)
      case( shell )
    end select
 
  end subroutine mat_c2d
 
  ! (Gaku Hashimoto, The University of Tokyo, 2012/11/15) <
  !####################################################################
  subroutine matlmatrix_shell                        &
      (gauss, secttype, d,                    &
      e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, &
      alpha, n_layer)
    !####################################################################
 
    type(tgaussstatus), intent(in)  :: gauss
    integer, intent(in)             :: sectType, n_layer
    real(kind = kreal), intent(out) :: d(:, :)
    real(kind = kreal), intent(in)  :: e1_hat(3), e2_hat(3), e3_hat(3)
    real(kind = kreal), intent(in)  :: cg1(3), cg2(3), cg3(3)
    real(kind = kreal), intent(out) :: alpha
 
    !--------------------------------------------------------------------
 
    real(kind = kreal)       :: c(3, 3, 3, 3)
    type(tmaterial), pointer :: matl
 
    !--------------------------------------------------------------------
 
    matl => gauss%pMaterial
 
    !--------------------------------------------------------------------
 
    if( iselastic(matl%mtype) ) then
      call linearelastic_shell                     &
        (matl, secttype, c,                     &
        e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, &
        alpha, n_layer)
 
      call mat_c2d_shell(c, d, secttype)
    else
      stop "Material type not supported!"
    end if
 
    !--------------------------------------------------------------------
 
    return
 
    !####################################################################
  end subroutine matlmatrix_shell
  !####################################################################
  ! > (Gaku Hashimoto, The University of Tokyo, 2012/11/15)
 
 
  ! (Gaku Hashimoto, The University of Tokyo, 2012/11/15) <
  !####################################################################
  subroutine mat_c2d_shell(c, D, itype)
    !####################################################################
 
    real(kind = kreal), intent(in)  :: c(:, :, :, :)
    real(kind = kreal), intent(out) :: d(:, :)
    integer, intent(in)             :: itype
 
    !--------------------------------------------------------------------
 
    integer :: index_i(5), index_j(5), &
      index_k(5), index_l(5)
    integer :: i, j, k, l
    integer :: is, js
 
    !--------------------------------------------------------------------
 
    index_i(1) = 1
    index_i(2) = 2
    index_i(3) = 1
    index_i(4) = 2
    index_i(5) = 3
 
    index_j(1) = 1
    index_j(2) = 2
    index_j(3) = 2
    index_j(4) = 3
    index_j(5) = 1
 
    index_k(1) = 1
    index_k(2) = 2
    index_k(3) = 1
    index_k(4) = 2
    index_k(5) = 3
 
    index_l(1) = 1
    index_l(2) = 2
    index_l(3) = 2
    index_l(4) = 3
    index_l(5) = 1
 
    !--------------------------------------------------------------------
 
    d(:, :) = 0.0d0
 
    !--------------------------------------------------------------------
 
    select case( itype )
      case( shell )
 
        do js = 1, 5
 
          do is = 1, 5
 
            i = index_i(is)
            j = index_j(is)
            k = index_k(js)
            l = index_l(js)
 
            d(is, js) = c(i, j, k, l)
 
          end do
 
        end do
 
    end select
 
    !--------------------------------------------------------------------
 
    return
 
    !####################################################################
  end subroutine mat_c2d_shell
  !####################################################################
  ! > (Gaku Hashimoto, The University of Tokyo, 2012/11/15)
 
  subroutine getconnectorproperty( gauss, ctype, dofid, rparams, iparams, stretch )
    type( tgaussstatus ), intent(in)         :: gauss
    integer(kind=kint), intent(in)           :: ctype
    integer(kind=kint), intent(in)           :: dofid
    real(kind=kreal), intent(out)            :: rparams(:)
    integer(kind=kint), intent(out)          :: iparams(:)
    real(kind=kreal), intent(in), optional   :: stretch
 
    real(kind=kreal) ina(1), outa(4)
    logical :: ierr
    character(len=DICT_KEY_LENGTH) :: cnkey
 
    rparams(:) = 0.d0
    iparams(:) = 0
 
    if( ctype == m_spring_dof ) then ! Dof Spring
      write(cnkey,'(A,I0)') trim(mc_spring),dofid
      iparams(1) = gauss%pMaterial%variables_i(m_spring_d_ndoffset+2*dofid  )
      iparams(2) = gauss%pMaterial%variables_i(m_spring_d_ndoffset+2*dofid+1)
    else if( ctype == m_spring_axial ) then ! Dof Spring
      write(cnkey,'(A,A)') trim(mc_spring),'_A'
    else if( ctype == m_dashpot_dof ) then ! Dof Dashpot
      write(cnkey,'(A,I0)') trim(mc_dashpot),dofid
      iparams(1) = gauss%pMaterial%variables_i(m_dashpot_d_ndoffset+2*dofid  )
      iparams(2) = gauss%pMaterial%variables_i(m_dashpot_d_ndoffset+2*dofid+1)
    else if( ctype == m_dashpot_axial ) then ! Dof Dashpot
      write(cnkey,'(A,A)') trim(mc_dashpot),'_A'
    else
      stop "CONNECTOR ctype is not defined"
    endif
 
    if( present(stretch) ) then
      ina(1) = stretch
      call fetch_tabledata( cnkey, gauss%pMaterial%dict, outa(1:1), ierr, ina )
    else
      call fetch_tabledata( cnkey, gauss%pMaterial%dict, outa(1:1), ierr )
    endif
    rparams(1) = outa(1)
 
  end subroutine
 
end module m_matmatrix