fstr_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_fstr_elemact
  use hecmw
  use m_fstr
  use m_elemact
 
contains
 
  subroutine apply_amplitude_control( hecMESH, elemact, dumid, elements, amp_id, ctime )
    type(hecmwst_local_mesh), intent(in)   :: hecMESH
    type(telemact), intent(in)             :: elemact
    integer(kind=kint), intent(in)         :: dumid, amp_id
    type(telement), pointer, intent(inout) :: elements(:)
    real(kind=kreal), intent(in)           :: ctime
    
    real(kind=kreal)   :: amp_val
    integer(kind=kint) :: amp_state
    
    call hecmw_get_amplitude_value(hecmesh%amp, amp_id, ctime, amp_val)
    
    ! Simple rule: amp_val > 0.5 -> ACTIVE, otherwise INACTIVE
    if( amp_val > 0.5d0 ) then
      amp_state = kelact_active
    else
      amp_state = kelact_inactive
    endif
    
    call set_elemact_flag( hecmesh, elemact, dumid, elements, amp_state, .false. )
  end subroutine apply_amplitude_control
 
  subroutine fstr_update_elemact_solid( hecMESH, fstrSOLID, cstep, ctime )
    type(hecmwst_local_mesh), intent(in)   :: hecMESH
    type(fstr_solid), intent(inout)        :: fstrSOLID
    integer(kind=kint), intent(in)         :: cstep
    real(kind=kreal), intent(in)           :: ctime        
 
    integer(kind=kint) :: idum, amp_id, gid
    real(kind=kreal)   :: amp_val
    integer(kind=kint) :: target_state  ! Target state from control file (STATE=ON/OFF)
 
    do idum = 1, fstrsolid%elemact%ELEMACT_egrp_tot
      gid = fstrsolid%elemact%ELEMACT_egrp_GRPID(idum)
      if( .not. fstr_iselemactivationactive( fstrsolid, gid, cstep ) ) then
        call set_elemact_flag( hecmesh, fstrsolid%elemact, idum, fstrsolid%elements, kelact_active, .false. )
        cycle
      endif
 
      amp_id = fstrsolid%elemact%ELEMACT_egrp_amp(idum)
      if( amp_id > 0 ) then
        ! Amplitude control overrides STATE
        call apply_amplitude_control( hecmesh, fstrsolid%elemact, idum, fstrsolid%elements, amp_id, ctime )
        cycle
      end if
 
      ! No amplitude: use STATE from control file
      target_state = fstrsolid%elemact%ELEMACT_egrp_state(idum)
 
      if( fstrsolid%elemact%ELEMACT_egrp_depends(idum) == kelactd_none ) then
        call set_elemact_flag( hecmesh, fstrsolid%elemact, idum, fstrsolid%elements, target_state, .false. )
      else
        call set_elemact_flag( hecmesh, fstrsolid%elemact, idum, fstrsolid%elements, target_state, .true. )
      endif
    end do
 
  end subroutine
 
  subroutine fstr_update_elemact_solid_by_value( hecMESH, fstrSOLID, cstep, ctime )
    type(hecmwst_local_mesh), intent(in)   :: hecMESH
    type(fstr_solid), intent(inout)        :: fstrSOLID
    integer(kind=kint), intent(in)         :: cstep
    real(kind=kreal), intent(in)           :: ctime        
 
    integer(kind=kint) :: idum, amp_id, gid
    integer(kind=kint) :: n_changed_local, n_changed_total
    real(kind=kreal)   :: amp_val
 
    n_changed_total = 0
 
    do idum = 1, fstrsolid%elemact%ELEMACT_egrp_tot
      gid = fstrsolid%elemact%ELEMACT_egrp_GRPID(idum)
      if( .not. fstr_iselemactivationactive( fstrsolid, gid, cstep ) ) cycle
 
      amp_id = fstrsolid%elemact%ELEMACT_egrp_amp(idum)
      if( amp_id > 0 ) then
        ! Amplitude control overrides stress-based control
        call apply_amplitude_control( hecmesh, fstrsolid%elemact, idum, fstrsolid%elements, amp_id, ctime )
        cycle
      end if
 
      ! No amplitude: use stress-based control
      n_changed_local = 0
      call activate_elemact_flag_by_value( hecmesh, fstrsolid%elemact, idum, fstrsolid%elements, n_changed_local )
      n_changed_total = n_changed_total + n_changed_local
    end do
 
    fstrsolid%elemact%ELEMACT_n_changed = n_changed_total
 
    if( hecmesh%my_rank == 0 .and. n_changed_total > 0 ) then
      write(*,'(a,i6,a)') ' *** ELEMACT: ', n_changed_total, ' element(s) changed state (ACTIVE -> INACTIVE)'
    endif
 
  end subroutine
 
  subroutine fstr_update_elemact_heat( hecMESH, elemact, ctime, elements )
    type(hecmwst_local_mesh), intent(in)   :: hecmesh
    type(telemact), intent(in)               :: elemact
    real(kind=kreal), intent(in)           :: ctime        
    type(telement), pointer, intent(inout) :: elements(:)
 
    integer(kind=kint) :: idum, amp_id
    real(kind=kreal)   :: amp_val
 
    call clear_elemact_flag_all( hecmesh, elemact, elements )
 
    do idum = 1, elemact%ELEMACT_egrp_tot
      amp_id = elemact%ELEMACT_egrp_amp(idum)
      amp_val = 1.d0
      if( amp_id > 0 ) then
        call hecmw_get_amplitude_value(hecmesh%amp, amp_id, ctime, amp_val)
        if( amp_val < 1.d0 ) cycle
      end if
 
      call activate_elemact_flag( hecmesh, elemact, idum, elements )
    end do
 
  end subroutine
 
  subroutine fstr_updatedof_elemact( ndof, hecMESH, elemact, elements, vec_old, vec_new )
    integer(kind=kint), intent(in)            :: ndof
    type(hecmwst_local_mesh), intent(in)      :: hecMESH
    type(telemact), intent(in)                  :: elemact
    type(telement), pointer, intent(inout)    :: elements(:)
    real(kind=kreal), pointer, intent(in)     :: vec_old(:)
    real(kind=kreal), pointer, intent(inout)  :: vec_new(:)
 
    integer(kind=kint) :: icel, in0
    integer(kind=kint) :: iS, iE, ic_type, nodlocal, i
    real(kind=kreal), pointer   :: active(:)
 
    allocate(active(hecmesh%n_node))
    active = -1.d0
 
    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(ic_type == 3414) cycle
 
      do icel = is, ie
        if( elements(icel)%elemact_flag == kelact_inactive ) cycle
        in0 = hecmesh%elem_node_index(icel-1)
        nn = hecmw_get_max_node(ic_type)
        do i = 1, nn
          nodlocal = hecmesh%elem_node_item(in0+i)
          active(nodlocal) = 1.d0
        enddo
      enddo
    enddo
 
    call hecmw_update_r(hecmesh,active,hecmesh%n_node,1)
 
    do i = 1, hecmesh%n_node
      if( active(i) > 0.d0 ) cycle
      vec_new(ndof*(i-1)+1:ndof*i) = vec_old(ndof*(i-1)+1:ndof*i)
    end do
 
    deallocate(active)
 
  end subroutine
 
  subroutine output_elemact_flag( hecMESH, elements, outval )
    type(hecmwst_local_mesh), intent(in)     :: hecmesh
    type(telement), pointer, intent(in)      :: elements(:)
    real(kind=kreal), pointer, intent(inout) :: outval(:)    
 
    integer(kind=kint) :: icel
 
    outval = 0.d0
    do icel = 1, hecmesh%n_elem
      if( elements(icel)%elemact_flag == kelact_inactive ) outval(icel) = 1.d0
    end do
 
  end subroutine
 
  subroutine activate_elemact_flag( hecMESH, elemact, dumid, elements )
    type(hecmwst_local_mesh), intent(in)   :: hecMESH
    type(telemact), intent(in)               :: elemact
    integer(kind=kint), intent(in)         :: dumid
    type(telement), pointer, intent(inout) :: elements(:)
 
    integer(kind=kint) :: ig, iS0, iE0, ik, icel
 
    if( dumid < 0 .or. dumid > elemact%ELEMACT_egrp_tot ) return
 
    ig = elemact%ELEMACT_egrp_ID(dumid)
    is0 = hecmesh%elem_group%grp_index(ig-1) + 1
    ie0 = hecmesh%elem_group%grp_index(ig  )
 
    do ik=is0,ie0
      icel = hecmesh%elem_group%grp_item(ik)
      elements(icel)%elemact_flag = kelact_inactive
      elements(icel)%elemact_coeff = elemact%ELEMACT_egrp_eps(dumid)
    end do
 
  end subroutine
 
  subroutine activate_elemact_flag_by_value( hecMESH, elemact, dumid, elements, n_changed )
    type(hecmwst_local_mesh), intent(in)   :: hecMESH
    type(telemact), intent(in)               :: elemact
    integer(kind=kint), intent(in)         :: dumid
    type(telement), pointer, intent(inout) :: elements(:)
    integer(kind=kint), intent(out)        :: n_changed
 
    integer(kind=kint) :: ig, iS0, iE0, ik, icel, dtype, ig0
    integer(kind=kint) :: old_flag
    real(kind=kreal)   :: thlow, thup, stress(6), mises, ps
    integer(kind=kint) :: target_state  ! Target state from control file (STATE=ON/OFF)
 
    n_changed = 0
    if( dumid < 0 .or. dumid > elemact%ELEMACT_egrp_tot ) return
    if( elemact%ELEMACT_egrp_depends(dumid) == kelactd_none ) return 
 
    ig = elemact%ELEMACT_egrp_ID(dumid)
    is0 = hecmesh%elem_group%grp_index(ig-1) + 1
    ie0 = hecmesh%elem_group%grp_index(ig  )
 
    thlow = elemact%ELEMACT_egrp_ts_lower(dumid)
    thup = elemact%ELEMACT_egrp_ts_upper(dumid)
 
    ! Get the target state from control file
    ! STATE=ON  -> kELACT_ACTIVE   (0): Elements start active
    ! STATE=OFF -> kELACT_INACTIVE (1): Elements start inactive
    target_state = kelact_inactive  ! Default is INACTIVE
    if (associated(elemact%ELEMACT_egrp_state)) then
      target_state = elemact%ELEMACT_egrp_state(dumid)
    endif
 
    do ik=is0,ie0
      icel = hecmesh%elem_group%grp_item(ik)
      
      ! Case 1: target_state is INACTIVE (STATE=OFF in control file)
      ! Set to INACTIVE regardless of stress (no stress dependency check needed)
      if( target_state == kelact_inactive ) then
        elements(icel)%elemact_flag = kelact_inactive
        elements(icel)%elemact_coeff = elemact%ELEMACT_egrp_eps(dumid)
        cycle
      endif
      
      ! Case 2: Element is already INACTIVE
      ! Once deactivated, element stays INACTIVE (no reactivation: one-way transition)
      if( elements(icel)%elemact_flag == kelact_inactive ) cycle
      
      ! Case 3: target_state is ACTIVE (STATE=ON) and element is currently ACTIVE
      ! Check stress to decide whether to deactivate: ACTIVE -> INACTIVE transition
      do ig0=1,size(elements(icel)%gausses)
        ! get mises
        if( elemact%ELEMACT_egrp_depends(dumid) == kelactd_stress ) then
          stress(1:6) = elements(icel)%gausses(ig0)%stress(1:6)
        elseif( elemact%ELEMACT_egrp_depends(dumid) == kelactd_strain ) then
          stress(1:6) = elements(icel)%gausses(ig0)%strain(1:6)
        else
          return
        endif
        ps = ( stress(1) + stress(2) + stress(3) ) / 3.0d0
        mises = 0.5d0 * ( (stress(1)-ps)**2 + (stress(2)-ps)**2 + (stress(3)-ps)**2 )
        mises = mises + stress(4)**2 + stress(5)**2 + stress(6)**2
        mises = dsqrt( 3.0d0 * mises )
        
        ! If stress is OUT OF RANGE [thlow, thup], deactivate the element
        old_flag = elements(icel)%elemact_flag
        if( .not. (thlow <= mises .and. mises <= thup) ) then
          elements(icel)%elemact_flag = kelact_inactive
          elements(icel)%elemact_coeff = elemact%ELEMACT_egrp_eps(dumid)
        endif
        ! Count state change
        if( elements(icel)%elemact_flag /= old_flag ) n_changed = n_changed + 1
        
        exit
      enddo
    end do
 
  end subroutine
 
  subroutine set_elemact_flag( hecMESH, elemact, dumid, elements, flag, init_only )
    type(hecmwst_local_mesh), intent(in)   :: hecMESH
    type(telemact), intent(in)               :: elemact
    integer(kind=kint), intent(in)         :: dumid
    type(telement), pointer, intent(inout) :: elements(:)
    integer(kind=kint), intent(in)         :: flag
    logical, intent(in)                    :: init_only
 
    integer(kind=kint) :: ig, iS0, iE0, ik, icel
 
    if( dumid < 0 .or. dumid > elemact%ELEMACT_egrp_tot ) return
 
    ig = elemact%ELEMACT_egrp_ID(dumid)
    is0 = hecmesh%elem_group%grp_index(ig-1) + 1
    ie0 = hecmesh%elem_group%grp_index(ig  )
 
    do ik=is0,ie0
      icel = hecmesh%elem_group%grp_item(ik)
      if( init_only ) then
        if( elements(icel)%elemact_flag == kelact_undefined ) elements(icel)%elemact_flag = flag
      else
        elements(icel)%elemact_flag = flag
      endif
      elements(icel)%elemact_coeff = elemact%ELEMACT_egrp_eps(dumid)
    end do
 
  end subroutine
 
  subroutine clear_elemact_flag_all( hecMESH, elemact, elements )
    type(hecmwst_local_mesh), intent(in)   :: hecMESH
    type(telemact), intent(in)               :: elemact
    type(telement), pointer, intent(inout) :: elements(:)
 
    integer(kind=kint) :: idum, ig, iS0, iE0, ik, icel
 
    do idum = 1, elemact%ELEMACT_egrp_tot
      ig = elemact%ELEMACT_egrp_GRPID(idum)
      is0 = hecmesh%elem_group%grp_index(ig-1) + 1
      ie0 = hecmesh%elem_group%grp_index(ig  )
 
      do ik=is0,ie0
        icel = hecmesh%elem_group%grp_item(ik)
        elements(icel)%elemact_flag = kelact_undefined
      end do
    end do
 
  end subroutine
 
end module m_fstr_elemact