m_element_activation.f90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
! Copyright (c) 2016 The University of Tokyo
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
module m_elemact
  use hecmw
  use mmechgauss
 
  implicit none
 
  type telemact
    integer(kind=kint) :: elemact_egrp_tot
    integer(kind=kint), pointer :: elemact_egrp_grpid  (:)  =>null()
    integer(kind=kint), pointer :: elemact_egrp_id     (:)  =>null()
    integer(kind=kint), pointer :: elemact_egrp_amp    (:)  =>null()
    real(kind=kreal), pointer   :: elemact_egrp_eps(:)  =>null()
    integer(kind=kint), pointer :: elemact_egrp_depends(:)  =>null()
    real(kind=kreal), pointer   :: elemact_egrp_ts_lower(:)  =>null()
    real(kind=kreal), pointer   :: elemact_egrp_ts_upper(:)  =>null()
    integer(kind=kint), pointer :: elemact_egrp_state  (:)  =>null()  ! Element activation state
    integer(kind=kint) :: elemact_n_changed = 0  ! Number of elements that changed state in last update
    logical :: elemact_conv_deferred = .false.   ! True if convergence was deferred due to state change (once per substep)
  end type
 
  ! Element status flags
  integer, parameter :: kelact_undefined = -1  ! Undefined state
  integer, parameter :: kelact_active  = 0   ! Element is active in calculation
  integer, parameter :: kelact_inactive = 1   ! Element is inactive (deactivated) in calculation
 
  ! Dependency type flags
  integer, parameter :: kelactd_none    = 1   ! No dependencies
  integer, parameter :: kelactd_stress  = 2   ! Stress dependent
  integer, parameter :: kelactd_strain  = 3   ! Strain dependent
 
contains
 
  
  subroutine print_elemact_info( elact )
    type(telemact), intent(in) :: elact
 
    integer(kind=kint) :: i
 
    write(*,'(A,I0)') 'ELEMACT_egrp_tot: ', elact%ELEMACT_egrp_tot
    do i=1,elact%ELEMACT_egrp_tot
      write(*,*) 'DUMMY : ',i
      write(*,*) 'GRPID    : ',elact%ELEMACT_egrp_GRPID
      write(*,*) 'EGRPID   : ',elact%ELEMACT_egrp_ID
      write(*,*) 'AMPLITUDE: ',elact%ELEMACT_egrp_amp
      write(*,*) 'EPSILON  : ',elact%ELEMACT_egrp_eps
      write(*,*) 'DEPENDS  : ',elact%ELEMACT_egrp_depends
      write(*,*) 'THLOWER  : ',elact%ELEMACT_egrp_ts_lower
      write(*,*) 'THUPPER  : ',elact%ELEMACT_egrp_ts_upper
    end do
  end subroutine
 
  
  subroutine stf_dummy( ndof, nn, ecoord, u, stiff, element )
    use mmechgauss
    integer(kind=kint), intent(in)  :: ndof
    integer(kind=kint), intent(in)  :: nn
    real(kind=kreal),   intent(in)  :: ecoord(3,nn)           
    real(kind=kreal), intent(in)    :: u(3,nn)                
    type(telement), intent(inout)   :: element
    real(kind=kreal),   intent(out) :: stiff(:,:)             
 
    integer(kind=kint) :: i,j,k,m,n
    real(kind=kreal)   :: dnn, coeff, xmax(3), xmin(3), dl
 
    !get average element size
    xmax(1:3) = ecoord(1:3,1)
    xmin(1:3) = ecoord(1:3,1)
    do i=2,nn
      do k=1,ndof
        if( ecoord(k,i) > xmax(k) ) xmax(k) = ecoord(k,i)
        if( ecoord(k,i) < xmin(k) ) xmin(k) = ecoord(k,i)
      end do
    end do
    dl = 0.33333333333d0*((xmax(1)-xmin(1))+(xmax(2)-xmin(2))+(xmax(3)-xmin(3)))
    if( dl < 1.d-8 ) dl = 1.d0
 
    coeff = element%elemact_coeff/dl
    dnn = 1.d0/dble(nn)
    stiff(:,:) = 0.d0
    do i=1,nn
      m = ndof*(i-1)
 
      do j=1,nn
        n = ndof*(j-1)
        do k=1,ndof
          stiff(m+k,n+k) = stiff(m+k,n+k)-coeff*dnn
        end do
      end do
 
      do k=1,ndof
        stiff(m+k,m+k) = stiff(m+k,m+k)+coeff
      end do
    end do
 
  end subroutine
 
  
  subroutine update_dummy( ndof, nn, ecoord, u, du, qf, element )
    use mmechgauss
    integer(kind=kint), intent(in)  :: ndof
    integer(kind=kint), intent(in)  :: nn
    real(kind=kreal),   intent(in)  :: ecoord(3,nn)           
    real(kind=kreal), intent(in)    :: u(3,nn)                
    real(kind=kreal), intent(in)    :: du(3,nn)               
    type(telement), intent(inout)   :: element
    real(kind=kreal),   intent(out) :: qf(:)                  
 
    integer(kind=kint) :: i,k,m
    real(kind=kreal)   :: dnn, coeff, aveu(3), xmax(3), xmin(3), dl
 
    !get average element size
    xmax(1:3) = ecoord(1:3,1)
    xmin(1:3) = ecoord(1:3,1)
    do i=2,nn
      do k=1,ndof
        if( ecoord(k,i) > xmax(k) ) xmax(k) = ecoord(k,i)
        if( ecoord(k,i) < xmin(k) ) xmin(k) = ecoord(k,i)
      end do
    end do
    dl = 0.33333333333d0*((xmax(1)-xmin(1))+(xmax(2)-xmin(2))+(xmax(3)-xmin(3)))
    if( dl < 1.d-8 ) dl = 1.d0
 
    coeff = element%elemact_coeff/dl
    dnn = 1.d0/dble(nn)
    aveu(:) = 0.d0
    do i=1,nn
      do k=1,ndof
        aveu(k) = aveu(k) + u(k,i) + du(k,i)
      end do
    end do
    aveu(:) = aveu(:)*dnn
 
    qf(:) = 0.d0
    do i=1,nn
      m = ndof*(i-1)
 
      do k=1,ndof
        qf(m+k) = coeff*(u(k,i)+du(k,i)-aveu(k))
      end do
    end do
 
    do i=1,size(element%gausses)
      element%gausses(i)%strain = 0.d0
      element%gausses(i)%stress = 0.d0
      element%gausses(i)%strain_out = 0.d0
      element%gausses(i)%stress_out = 0.d0
      if( associated(element%gausses(i)%istatus) ) element%gausses(i)%istatus = 0
      if( associated(element%gausses(i)%fstatus) ) element%gausses(i)%fstatus = 0.d0
    end do
 
  end subroutine
 
 
end module m_elemact