fstr_NodalKinematics.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_nodalkinematics
  use m_fstr
  use m_fstr_finiterotationkinematics, only: fstr_uses_finite_rotation_kinematics, &
    shellorthonormalizetriad, shellupdatetriadwithincrement, shellcomposerotationvector
  implicit none
 
  private
 
  public :: fstr_ensure_finite_rotation_state
  public :: fstr_begin_nodal_kinematics_step
  public :: fstr_apply_solution_increment
  public :: fstr_commit_solution_increment
 
contains
 
  subroutine fstr_ensure_finite_rotation_state( hecMESH, fstrSOLID, ndof )
    use elementinfo, only: fe_mitc4_shell
    implicit none
 
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (fstr_solid), intent(inout)      :: fstrsolid
    integer(kind=kint), intent(in)        :: ndof
 
    integer(kind=kint) :: itype, is, ie, ic_type, icel, iis, nn, j, node_id
    integer(kind=kint), allocatable :: node_mode(:), node_count(:)
    real(kind=kreal), allocatable :: director_sum(:,:), tangent_sum(:,:)
    real(kind=kreal) :: ecoord(3, 8), triad(3, 3), trial(3, 3)
    real(kind=kreal) :: director(3), tangent(3), ref_axis(3), normv, proj
 
    if( .not. fstrsolid%has_finite_rotation_kinematics ) return
    if( fstrsolid%finite_rotation_state_ready ) return
    if( .not. associated(fstrsolid%shell_rot_state) ) return
 
    allocate( director_sum(3, hecmesh%n_node) )
    allocate( tangent_sum(3, hecmesh%n_node) )
    allocate( node_mode(hecmesh%n_node) )
    allocate( node_count(hecmesh%n_node) )
    director_sum(:, :) = 0.0d0
    tangent_sum(:, :) = 0.0d0
    node_mode(:) = 0
    node_count(:) = 0
 
    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 /= fe_mitc4_shell ) cycle
      do icel = is, ie
        iis = hecmesh%elem_node_index(icel-1)
        nn = hecmesh%elem_node_index(icel) - iis
        if( .not. associated( fstrsolid%elements(icel)%gausses ) ) cycle
        if( .not. fstr_uses_finite_rotation_kinematics( ic_type, nn, &
            fstrsolid%elements(icel)%gausses(1)%pMaterial ) ) cycle
 
        do j = 1, min(nn, 8)
          node_id = hecmesh%elem_node_item(iis+j)
          ecoord(1:3, j) = hecmesh%node(3*node_id-2:3*node_id)
        end do
 
        if( ndof >= 6 ) then
          do j = 1, nn
            call fstr_reference_shell_triad( nn, ecoord(1:3, 1:nn), j, triad )
            node_id = hecmesh%elem_node_item(iis+j)
            if( node_id <= 0 .or. node_id > hecmesh%n_node ) cycle
            if( fstrsolid%shell_rot_state(node_id) /= 0 ) cycle
            ! average shell triads at shared nodes
            director_sum(1:3, node_id) = director_sum(1:3, node_id) + triad(1:3, 3)
            tangent_sum(1:3, node_id) = tangent_sum(1:3, node_id) + triad(1:3, 1)
            node_count(node_id) = node_count(node_id) + 1
            node_mode(node_id) = 1
          end do
        endif
      end do
    end do
 
    do node_id = 1, hecmesh%n_node
      if( node_count(node_id) <= 0 ) cycle
      if( fstrsolid%shell_rot_state(node_id) /= 0 ) cycle
 
      director(1:3) = director_sum(1:3, node_id)
      normv = dsqrt( dot_product( director(1:3), director(1:3) ) )
      if( normv <= 1.0d-14 ) then
        call fstr_set_identity_triad( triad )
        call fstr_store_shell_triad_node( fstrsolid, node_id, triad, node_mode(node_id) )
        cycle
      endif
      director(1:3) = director(1:3)/normv
 
      tangent(1:3) = tangent_sum(1:3, node_id)
      proj = dot_product( tangent(1:3), director(1:3) )
      tangent(1:3) = tangent(1:3) - proj*director(1:3)
      normv = dsqrt( dot_product( tangent(1:3), tangent(1:3) ) )
      if( normv <= 1.0d-14 ) then
        ref_axis(1:3) = (/ 1.0d0, 0.0d0, 0.0d0 /)
        if( dabs(director(1)) > 0.9d0 ) ref_axis(1:3) = (/ 0.0d0, 1.0d0, 0.0d0 /)
        proj = dot_product( ref_axis(1:3), director(1:3) )
        tangent(1:3) = ref_axis(1:3) - proj*director(1:3)
        normv = dsqrt( dot_product( tangent(1:3), tangent(1:3) ) )
      endif
      tangent(1:3) = tangent(1:3)/normv
 
      call fstr_set_identity_triad( trial )
      trial(1:3, 1) = tangent(1:3)
      trial(1:3, 3) = director(1:3)
      call shellorthonormalizetriad( trial, triad )
      call fstr_store_shell_triad_node( fstrsolid, node_id, triad, node_mode(node_id) )
    end do
 
    deallocate( director_sum )
    deallocate( tangent_sum )
    deallocate( node_mode )
    deallocate( node_count )
    fstrsolid%finite_rotation_state_ready = .true.
 
  end subroutine fstr_ensure_finite_rotation_state
 
  subroutine fstr_begin_nodal_kinematics_step( hecMESH, fstrSOLID, ndof )
    implicit none
 
    type (hecmwst_local_mesh), intent(in) :: hecmesh
    type (fstr_solid), intent(inout)      :: fstrsolid
    integer(kind=kint), intent(in)        :: ndof
 
    call fstr_ensure_finite_rotation_state( hecmesh, fstrsolid, ndof )
    if( .not. associated(fstrsolid%shell_triad) ) return
 
    fstrsolid%shell_triad_bak(:) = fstrsolid%shell_triad(:)
    fstrsolid%shell_drill_bak(:) = fstrsolid%shell_drill(:)
    fstrsolid%shell_dtriad(:) = fstrsolid%shell_triad(:)
    fstrsolid%shell_ddrill(:) = fstrsolid%shell_drill(:)
 
  end subroutine fstr_begin_nodal_kinematics_step
 
  subroutine fstr_apply_solution_increment( hecMESH, fstrSOLID, ndof, x )
    implicit none
 
    type (hecmwst_local_mesh), intent(in)    :: hecmesh
    type (fstr_solid), intent(inout)         :: fstrsolid
    integer(kind=kint), intent(in)           :: ndof
    real(kind=kreal), intent(in)             :: x(:)
 
    integer(kind=kint) :: node_id, idx, base
    real(kind=kreal) :: theta_inc(3), theta_compat(3)
    real(kind=kreal) :: triad_old(3, 3), triad_new(3, 3), drill_new
 
    if( .not. fstrsolid%has_finite_rotation_kinematics ) then
      do node_id = 1, hecmesh%n_node
        idx = ndof*(node_id-1)
        fstrsolid%dunode(idx+1:idx+ndof) = fstrsolid%dunode(idx+1:idx+ndof) + x(idx+1:idx+ndof)
      end do
      return
    endif
    call fstr_ensure_finite_rotation_state( hecmesh, fstrsolid, ndof )
    if( .not. associated(fstrsolid%shell_node_mode) ) then
      do node_id = 1, hecmesh%n_node
        idx = ndof*(node_id-1)
        fstrsolid%dunode(idx+1:idx+ndof) = fstrsolid%dunode(idx+1:idx+ndof) + x(idx+1:idx+ndof)
      end do
      return
    endif
 
    do node_id = 1, hecmesh%n_node
      idx = ndof*(node_id-1)
      if( fstrsolid%shell_node_mode(node_id) == 1 ) then
        fstrsolid%dunode(idx+1:idx+3) = fstrsolid%dunode(idx+1:idx+3) + x(idx+1:idx+3)
        theta_inc(1:3) = x(idx+4:idx+6)
        base = 9*(node_id-1)
        triad_old(1:3, 1) = fstrsolid%shell_dtriad(base+1:base+3)
        triad_old(1:3, 2) = fstrsolid%shell_dtriad(base+4:base+6)
        triad_old(1:3, 3) = fstrsolid%shell_dtriad(base+7:base+9)
        call shellupdatetriadwithincrement( triad_old, fstrsolid%shell_ddrill(node_id), &
          theta_inc, triad_new, drill_new )
        fstrsolid%shell_dtriad(base+1:base+3) = triad_new(1:3, 1)
        fstrsolid%shell_dtriad(base+4:base+6) = triad_new(1:3, 2)
        fstrsolid%shell_dtriad(base+7:base+9) = triad_new(1:3, 3)
        fstrsolid%shell_ddrill(node_id) = drill_new
        ! update nodal rotation vector for output
        call shellcomposerotationvector( fstrsolid%dunode(idx+4:idx+6), x(idx+4:idx+6), theta_compat )
        fstrsolid%dunode(idx+4:idx+6) = theta_compat(1:3)
        if( ndof > 6 ) then
          fstrsolid%dunode(idx+7:idx+ndof) = fstrsolid%dunode(idx+7:idx+ndof) + x(idx+7:idx+ndof)
        endif
      else
        fstrsolid%dunode(idx+1:idx+ndof) = fstrsolid%dunode(idx+1:idx+ndof) + x(idx+1:idx+ndof)
      endif
    end do
 
  end subroutine fstr_apply_solution_increment
 
  subroutine fstr_commit_solution_increment( hecMESH, fstrSOLID, ndof )
    implicit none
 
    type (hecmwst_local_mesh), intent(in)    :: hecmesh
    type (fstr_solid), intent(inout)         :: fstrsolid
    integer(kind=kint), intent(in)           :: ndof
 
    integer(kind=kint) :: node_id, idx, base
    real(kind=kreal) :: theta_compat(3)
 
    if( .not. fstrsolid%has_finite_rotation_kinematics ) then
      do node_id = 1, hecmesh%n_node
        idx = ndof*(node_id-1)
        fstrsolid%unode(idx+1:idx+ndof) = fstrsolid%unode(idx+1:idx+ndof) + fstrsolid%dunode(idx+1:idx+ndof)
      end do
      return
    endif
    call fstr_ensure_finite_rotation_state( hecmesh, fstrsolid, ndof )
    if( .not. associated(fstrsolid%shell_node_mode) ) then
      do node_id = 1, hecmesh%n_node
        idx = ndof*(node_id-1)
        fstrsolid%unode(idx+1:idx+ndof) = fstrsolid%unode(idx+1:idx+ndof) + fstrsolid%dunode(idx+1:idx+ndof)
      end do
      return
    endif
 
    do node_id = 1, hecmesh%n_node
      idx = ndof*(node_id-1)
      if( fstrsolid%shell_node_mode(node_id) == 1 ) then
        fstrsolid%unode(idx+1:idx+3) = fstrsolid%unode(idx+1:idx+3) + fstrsolid%dunode(idx+1:idx+3)
        base = 9*(node_id-1)
        fstrsolid%shell_triad(base+1:base+9) = fstrsolid%shell_dtriad(base+1:base+9)
        fstrsolid%shell_drill(node_id) = fstrsolid%shell_ddrill(node_id)
        ! update nodal rotation vector for output
        call shellcomposerotationvector( fstrsolid%unode(idx+4:idx+6), fstrsolid%dunode(idx+4:idx+6), theta_compat )
        fstrsolid%unode(idx+4:idx+6) = theta_compat(1:3)
        if( ndof > 6 ) then
          fstrsolid%unode(idx+7:idx+ndof) = fstrsolid%unode(idx+7:idx+ndof) + fstrsolid%dunode(idx+7:idx+ndof)
        endif
      else
        fstrsolid%unode(idx+1:idx+ndof) = fstrsolid%unode(idx+1:idx+ndof) + fstrsolid%dunode(idx+1:idx+ndof)
      endif
    end do
 
  end subroutine fstr_commit_solution_increment
 
  ! ---------------------------------------------------------------------------
  ! Private helpers: reference-frame construction at element nodes.
  ! ---------------------------------------------------------------------------
 
  subroutine fstr_set_identity_triad( triad )
    implicit none
 
    real(kind=kreal), intent(out) :: triad(3, 3)
 
    triad(:, :) = 0.0d0
    triad(1, 1) = 1.0d0
    triad(2, 2) = 1.0d0
    triad(3, 3) = 1.0d0
 
  end subroutine fstr_set_identity_triad
 
  subroutine fstr_reference_shell_triad( nn, ecoord, inode, triad )
    use elementinfo, only: fe_mitc4_shell, getshapederiv
    implicit none
 
    integer(kind=kint), intent(in) :: nn
    integer(kind=kint), intent(in) :: inode
    real(kind=kreal), intent(in)  :: ecoord(3, nn)
    real(kind=kreal), intent(out) :: triad(3, 3)
 
    real(kind=kreal) :: xi, eta
    real(kind=kreal) :: shapederiv(4, 2)
    real(kind=kreal) :: g1(3), g2(3), e0(3), trial(3, 3), normv
    integer(kind=kint) :: i
 
    call fstr_set_identity_triad( trial )
    if( nn /= 4 ) then
      triad(1:3, 1:3) = trial(1:3, 1:3)
      return
    endif
 
    call getshapederiv( fe_mitc4_shell, (/ 0.0d0, 0.0d0 /), shapederiv )
    e0(1:3) = 0.0d0
    do i = 1, 4
      e0(1:3) = e0(1:3) + shapederiv(i, 1)*ecoord(1:3, i)
    end do
 
    select case( inode )
    case( 1 )
      xi = -1.0d0
      eta = -1.0d0
    case( 2 )
      xi =  1.0d0
      eta = -1.0d0
    case( 3 )
      xi =  1.0d0
      eta =  1.0d0
    case default
      xi = -1.0d0
      eta =  1.0d0
    end select
 
    call getshapederiv( fe_mitc4_shell, (/ xi, eta /), shapederiv )
    g1(1:3) = 0.0d0
    g2(1:3) = 0.0d0
    do i = 1, 4
      g1(1:3) = g1(1:3) + shapederiv(i, 1)*ecoord(1:3, i)
      g2(1:3) = g2(1:3) + shapederiv(i, 2)*ecoord(1:3, i)
    end do
 
    trial(1, 3) = g1(2)*g2(3) - g1(3)*g2(2)
    trial(2, 3) = g1(3)*g2(1) - g1(1)*g2(3)
    trial(3, 3) = g1(1)*g2(2) - g1(2)*g2(1)
 
    trial(1, 2) = trial(2, 3)*e0(3) - trial(3, 3)*e0(2)
    trial(2, 2) = trial(3, 3)*e0(1) - trial(1, 3)*e0(3)
    trial(3, 2) = trial(1, 3)*e0(2) - trial(2, 3)*e0(1)
    normv = dsqrt( dot_product( trial(1:3, 2), trial(1:3, 2) ) )
    if( normv > 1.0d-14 ) trial(1:3, 2) = trial(1:3, 2)/normv
 
    trial(1, 1) = trial(2, 2)*trial(3, 3) - trial(3, 2)*trial(2, 3)
    trial(2, 1) = trial(3, 2)*trial(1, 3) - trial(1, 2)*trial(3, 3)
    trial(3, 1) = trial(1, 2)*trial(2, 3) - trial(2, 2)*trial(1, 3)
 
    call shellorthonormalizetriad( trial, triad )
 
  end subroutine fstr_reference_shell_triad
 
  subroutine fstr_store_shell_triad_node( fstrSOLID, node_id, triad, mode )
    implicit none
 
    type (fstr_solid), intent(inout)  :: fstrSOLID
    integer(kind=kint), intent(in)    :: node_id, mode
    real(kind=kreal), intent(in)      :: triad(3, 3)
 
    integer(kind=kint) :: base
 
    if( node_id <= 0 ) return
    if( node_id > size(fstrsolid%shell_rot_state) ) return
    if( fstrsolid%shell_rot_state(node_id) /= 0 ) return
 
    base = 9*(node_id-1)
    fstrsolid%shell_rot_state(node_id) = mode
    ! initialize reference and current shell triads
    fstrsolid%shell_ref_triad(base+1:base+3) = triad(1:3, 1)
    fstrsolid%shell_ref_triad(base+4:base+6) = triad(1:3, 2)
    fstrsolid%shell_ref_triad(base+7:base+9) = triad(1:3, 3)
    fstrsolid%shell_triad(base+1:base+3) = triad(1:3, 1)
    fstrsolid%shell_triad(base+4:base+6) = triad(1:3, 2)
    fstrsolid%shell_triad(base+7:base+9) = triad(1:3, 3)
    fstrsolid%shell_triad_bak(base+1:base+9) = fstrsolid%shell_triad(base+1:base+9)
    fstrsolid%shell_dtriad(base+1:base+9) = fstrsolid%shell_triad(base+1:base+9)
    fstrsolid%shell_drill(node_id) = 0.0d0
    fstrsolid%shell_drill_bak(node_id) = 0.0d0
    fstrsolid%shell_ddrill(node_id) = 0.0d0
 
  end subroutine fstr_store_shell_triad_node
 
end module m_fstr_nodalkinematics