static_LIB_shell.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_static_lib_shell
  use hecmw, only : kint, kreal, hecmw_abort, hecmw_comm_get_comm
  use elementinfo
  use mitc_tying, only: numoftyingsets, numoftyingpoints, gettyingpoint, mitc9_xi_sign, mitc9_eta_sign
  use m_utilities, only: cross_product, calinverse
  use m_fstr_finiterotationkinematics, only: fstr_is_finite_rotation_shell_element, &
    shellrotationvectortomatrix, shellcomposerotationvector, shellskewmatrix
  use mmechgauss, only: fstr_uses_plane_stress_j2_shell
 
  implicit none
 
  private
 
  integer(kind=kint), parameter :: MAX_TYING_TERMS = 64
  integer(kind=kint), parameter :: SHELL_XI = 1_kint
  integer(kind=kint), parameter :: SHELL_ETA = 2_kint
  integer(kind=kint), parameter :: SHELL_ZETA = 3_kint
 
 
  public :: stf_shell_mitc
  public :: stf_shell_mitc33
  public :: elementstress_shell_mitc
  public :: dl_shell
  public :: dl_shell_33
  public :: update_shell_mitc
  public :: update_shell_mitc33
  public :: mass_shell
 
  !--------------------------------------------------------------------
 
  !--------------------------------------------------------------
  !
  ! (Programmer)
  ! Gaku Hashimoto
  ! Department of Human and Engineered Environmental Studies
  ! Graduate School of Frontier Sciences, The University of Tokyo
  ! 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8563 JAPAN
  !
  ! (Ref.)
  ! [1] Noguchi, H. and Hisada, T., "Sensitivity analysis in post-buckling problems of shell structures,"
  !     Computers & Structures, Vol.47, No.4, pp.699-710, (1993).
  ! [2] Dvorkin, E.N. and Bathe, K.J., "A Continuum Mechanics Based Four-node Shell Element for General Non-linear Analysis,"
  !     Engineering Computations, Vol.1, pp.77-88, (1984).
  ! [3] Bucalem, M.L. and Bathe, K.J., "Higher-order MITC general shell element,"
  !     International Journal for Numerical Methods in Engineering, Vol.36, pp.3729-3754, (1993).
  ! [4] Lee, P.S. and Bathe, K.J., "Development of MITC Isotropic Triangular Shell Finite Elements,"
  !     Computers & Structures, Vol.82, pp.945-962, (2004).
  !
  ! Xi YUAN
  !   Apr. 13, 2019: Introduce mass matrix calculation
  ! (Ref.)
  ! [5] E. Hinton, T. A. Rock and O. C. Zienkiewicz, "A Note on Mass Lumping and Related Processes in FEM,"
  !     Earthquake Engineering & Structural Dynamics, Vol.4, pp.245-249, (1976).
  !
  !--------------------------------------------------------------
 
contains
 
  pure function outer_product3(a, b) result(ab)
 
    real(kind = kreal), intent(in) :: a(3), b(3)
    real(kind = kreal) :: ab(3, 3)
    integer :: i, j
 
    do j = 1, 3
      do i = 1, 3
        ab(i, j) = a(i)*b(j)
      end do
    end do
 
  end function outer_product3
 
  !--------------------------------------------------------------------
  pure function shellmitc_symmetriccomponents(matrix) result(components)
    real(kind=kreal), intent(in) :: matrix(3, 3)
    real(kind=kreal) :: components(5)
 
    components = (/ matrix(1, 1), matrix(2, 2), matrix(1, 2)+matrix(2, 1), matrix(2, 3)+matrix(3, 2), &
      matrix(3, 1)+matrix(1, 3) /)
  end function shellmitc_symmetriccomponents
 
  !--------------------------------------------------------------------
  pure real(kind=kreal) function shellmitc_covariantjacobian(covariant_basis)
    real(kind=kreal), intent(in) :: covariant_basis(3, 3)
 
    shellmitc_covariantjacobian = covariant_basis(1, shell_xi) &
      *(covariant_basis(2, shell_eta)*covariant_basis(3, shell_zeta) &
      -covariant_basis(3, shell_eta)*covariant_basis(2, shell_zeta)) +covariant_basis(2, shell_xi) &
      *(covariant_basis(3, shell_eta)*covariant_basis(1, shell_zeta) &
      -covariant_basis(1, shell_eta)*covariant_basis(3, shell_zeta)) +covariant_basis(3, shell_xi) &
      *(covariant_basis(1, shell_eta)*covariant_basis(2, shell_zeta) &
      -covariant_basis(2, shell_eta)*covariant_basis(1, shell_zeta))
  end function shellmitc_covariantjacobian
 
  pure real(kind=kreal) function shellmitc_tyingzeta(etype, zeta)
    integer(kind=kint), intent(in) :: etype
    real(kind=kreal), intent(in) :: zeta
 
    shellmitc_tyingzeta = 0.0d0
    if( etype == fe_mitc9_shell ) shellmitc_tyingzeta = zeta
  end function shellmitc_tyingzeta
 
  !--------------------------------------------------------------------
  subroutine shellmitc_resolveformulation(etype, nn, ndof, gauss, has_nodal_state, &
      has_element_state, require_layer_state, kinematics, ndof_shell, finite_rotation, &
      use_director_tangent, use_green_lagrange, add_geometric_stiffness, update_state)
    use mmechgauss, only: tgaussstatus
    use mmaterial, only: infinitesimal, totallag, updatelag, iselastic, iselastoplastic
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn, ndof
    type(tGaussStatus), intent(in) :: gauss
    logical, intent(in) :: has_nodal_state, has_element_state, require_layer_state
    integer(kind=kint), intent(out) :: kinematics, ndof_shell
    logical, intent(out) :: finite_rotation, use_director_tangent, use_green_lagrange
    logical, intent(out) :: add_geometric_stiffness, update_state
    logical :: plane_stress_j2
 
    kinematics = gauss%pMaterial%nlgeom_flag
    if( .not. has_nodal_state ) kinematics = infinitesimal
    ndof_shell = min(ndof, 6_kint)
    finite_rotation = (kinematics == totallag .or. kinematics == updatelag) &
      .and. ndof_shell >= 6 .and. fstr_is_finite_rotation_shell_element(etype, nn) .and. iselastic(gauss%pMaterial%mtype)
    plane_stress_j2 = fstr_uses_plane_stress_j2_shell(etype, nn, gauss%pMaterial)
    use_director_tangent = finite_rotation
    use_green_lagrange = kinematics == totallag .and. finite_rotation
    add_geometric_stiffness = kinematics /= infinitesimal
    update_state = (finite_rotation .or. plane_stress_j2) .and. has_element_state
 
    if( kinematics /= infinitesimal ) then
      if( .not. finite_rotation ) call shellmitc_abortnonlinearunsupported(etype)
      if( require_layer_state .and. .not. update_state ) then
        call shellmitc_abortnonlinearunsupported(etype)
      endif
    endif
    if( iselastoplastic(gauss%pMaterial%mtype) .and. .not.plane_stress_j2 ) then
      call shellmitc_abortelastoplasticunsupported(etype)
    endif
  end subroutine shellmitc_resolveformulation
 
  !--------------------------------------------------------------------
  subroutine shellmitc_preparenodalkinematics(etype, nn, thick, kinematics, finite_rotation, &
      use_director_tangent, need_second_tangent, ecoord, nodal_state, ndtriad, ndreftriad, &
      ndcurtriad, evaluation_coords, v1, v2, v3, director, reference_director, director_tangent, director_second_tangent)
    use mmaterial, only: updatelag
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn, kinematics
    real(kind=kreal), intent(in) :: thick, ecoord(3, nn), nodal_state(6, nn)
    logical, intent(in) :: finite_rotation, use_director_tangent, need_second_tangent
    real(kind=kreal), intent(in), optional :: ndtriad(9, nn), ndreftriad(9, nn), ndcurtriad(9, nn)
    real(kind=kreal), intent(out) :: evaluation_coords(3, nn)
    real(kind=kreal), intent(out) :: v1(3, nn), v2(3, nn), v3(3, nn)
    real(kind=kreal), intent(out) :: director(3, nn), reference_director(3, nn)
    real(kind=kreal), intent(out) :: director_tangent(3, 3, nn)
    real(kind=kreal), intent(out) :: director_second_tangent(3, 3, 3, nn)
 
    evaluation_coords = ecoord
    if( kinematics == updatelag ) evaluation_coords = evaluation_coords+nodal_state(1:3, :)
    call shellmitc_setupnodaldirectors(etype, nn, thick, kinematics, evaluation_coords, &
      nodal_state, finite_rotation, use_director_tangent, need_second_tangent, ndtriad, &
      ndreftriad, ndcurtriad, v1, v2, v3, director, reference_director, director_tangent, director_second_tangent)
  end subroutine shellmitc_preparenodalkinematics
 
  !--------------------------------------------------------------------
  subroutine shellmitc_covariantbasis(nn, coords, director, zeta, shapefunc, shapederiv, covariant_basis)
    implicit none
 
    integer(kind=kint), intent(in) :: nn
    real(kind=kreal), intent(in) :: coords(3, nn), director(3, nn)
    real(kind=kreal), intent(in) :: zeta, shapefunc(nn), shapederiv(nn, 2)
    real(kind=kreal), intent(out) :: covariant_basis(3, 3)
 
    covariant_basis(:, shell_xi:shell_eta) = matmul(coords+zeta*director, shapederiv)
    covariant_basis(:, shell_zeta) = matmul(director, shapefunc)
 
  end subroutine shellmitc_covariantbasis
 
  !--------------------------------------------------------------------
  subroutine shellmitc_directorcontributions(nn, director, zeta, shapefunc, shapederiv, director_contribution)
    implicit none
 
    integer(kind=kint), intent(in) :: nn
    real(kind=kreal), intent(in) :: director(3, nn)
    real(kind=kreal), intent(in) :: zeta, shapefunc(nn), shapederiv(nn, 2)
    real(kind=kreal), intent(out) :: director_contribution(3, nn, 3)
 
    integer :: node
 
    do node = 1, nn
      director_contribution(:, node, shell_xi) = zeta*shapederiv(node, shell_xi)*director(:, node)
      director_contribution(:, node, shell_eta) = zeta*shapederiv(node, shell_eta)*director(:, node)
      director_contribution(:, node, shell_zeta) = shapefunc(node)*director(:, node)
    end do
  end subroutine shellmitc_directorcontributions
  pure subroutine shellmitc_basisfromcovariant(covariant_basis, local_basis, reciprocal_basis, det)
    real(kind=kreal), intent(in) :: covariant_basis(3, 3)
    real(kind=kreal), intent(out) :: local_basis(3, 3), reciprocal_basis(3, 3)
    real(kind=kreal), intent(out) :: det
 
    real(kind=kreal) :: det_inv, basis_norm
 
    det = shellmitc_covariantjacobian(covariant_basis)
    det_inv = 1.0d0/det
 
    reciprocal_basis(1, shell_xi) = det_inv *(covariant_basis(2, shell_eta)*covariant_basis(3, shell_zeta) &
      -covariant_basis(3, shell_eta)*covariant_basis(2, shell_zeta))
    reciprocal_basis(2, shell_xi) = det_inv *(covariant_basis(3, shell_eta)*covariant_basis(1, shell_zeta) &
      -covariant_basis(1, shell_eta)*covariant_basis(3, shell_zeta))
    reciprocal_basis(3, shell_xi) = det_inv *(covariant_basis(1, shell_eta)*covariant_basis(2, shell_zeta) &
      -covariant_basis(2, shell_eta)*covariant_basis(1, shell_zeta))
    reciprocal_basis(1, shell_eta) = det_inv *(covariant_basis(2, shell_zeta)*covariant_basis(3, shell_xi) &
      -covariant_basis(3, shell_zeta)*covariant_basis(2, shell_xi))
    reciprocal_basis(2, shell_eta) = det_inv *(covariant_basis(3, shell_zeta)*covariant_basis(1, shell_xi) &
      -covariant_basis(1, shell_zeta)*covariant_basis(3, shell_xi))
    reciprocal_basis(3, shell_eta) = det_inv *(covariant_basis(1, shell_zeta)*covariant_basis(2, shell_xi) &
      -covariant_basis(2, shell_zeta)*covariant_basis(1, shell_xi))
    reciprocal_basis(1, shell_zeta) = det_inv *(covariant_basis(2, shell_xi)*covariant_basis(3, shell_eta) &
      -covariant_basis(3, shell_xi)*covariant_basis(2, shell_eta))
    reciprocal_basis(2, shell_zeta) = det_inv *(covariant_basis(3, shell_xi)*covariant_basis(1, shell_eta) &
      -covariant_basis(1, shell_xi)*covariant_basis(3, shell_eta))
    reciprocal_basis(3, shell_zeta) = det_inv *(covariant_basis(1, shell_xi)*covariant_basis(2, shell_eta) &
      -covariant_basis(2, shell_xi)*covariant_basis(1, shell_eta))
 
    basis_norm = dsqrt(dot_product(covariant_basis(:, shell_zeta), covariant_basis(:, shell_zeta)))
    local_basis(:, shell_zeta) = covariant_basis(:, shell_zeta)/basis_norm
    local_basis(1, shell_xi) = covariant_basis(2, shell_eta)*local_basis(3, shell_zeta) &
      -covariant_basis(3, shell_eta)*local_basis(2, shell_zeta)
    local_basis(2, shell_xi) = covariant_basis(3, shell_eta)*local_basis(1, shell_zeta) &
      -covariant_basis(1, shell_eta)*local_basis(3, shell_zeta)
    local_basis(3, shell_xi) = covariant_basis(1, shell_eta)*local_basis(2, shell_zeta) &
      -covariant_basis(2, shell_eta)*local_basis(1, shell_zeta)
    basis_norm = dsqrt(dot_product(local_basis(:, shell_xi), local_basis(:, shell_xi)))
    local_basis(:, shell_xi) = local_basis(:, shell_xi)/basis_norm
    local_basis(1, shell_eta) = local_basis(2, shell_zeta)*local_basis(3, shell_xi) &
      -local_basis(3, shell_zeta)*local_basis(2, shell_xi)
    local_basis(2, shell_eta) = local_basis(3, shell_zeta)*local_basis(1, shell_xi) &
      -local_basis(1, shell_zeta)*local_basis(3, shell_xi)
    local_basis(3, shell_eta) = local_basis(1, shell_zeta)*local_basis(2, shell_xi) &
      -local_basis(2, shell_zeta)*local_basis(1, shell_xi)
    basis_norm = dsqrt(dot_product(local_basis(:, shell_eta), local_basis(:, shell_eta)))
    local_basis(:, shell_eta) = local_basis(:, shell_eta)/basis_norm
  end subroutine shellmitc_basisfromcovariant
 
  !--------------------------------------------------------------------
  subroutine shellmitc_preparepointkinematics(etype, nn, use_green_lagrange, reference_coords, &
      evaluation_coords, translation, director, reference_director, zeta, shapefunc, &
      shapederiv, covariant_basis, tangent_basis, reciprocal_basis, local_basis, &
      material_reciprocal_basis, material_local_basis, integration_jacobian, director_contribution)
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn
    logical, intent(in) :: use_green_lagrange
    real(kind=kreal), intent(in) :: reference_coords(3, nn), evaluation_coords(3, nn)
    real(kind=kreal), intent(in) :: translation(3, nn), director(3, nn)
    real(kind=kreal), intent(in) :: reference_director(3, nn)
    real(kind=kreal), intent(in) :: zeta, shapefunc(nn), shapederiv(nn, 2)
    real(kind=kreal), intent(out) :: covariant_basis(3, 3), tangent_basis(3, 2)
    real(kind=kreal), intent(out) :: reciprocal_basis(3, 3), local_basis(3, 3)
    real(kind=kreal), intent(out) :: material_reciprocal_basis(3, 3)
    real(kind=kreal), intent(out) :: material_local_basis(3, 3)
    real(kind=kreal), intent(out) :: integration_jacobian
    real(kind=kreal), intent(out) :: director_contribution(3, nn, 3)
 
    real(kind=kreal) :: reference_basis(3, 3), translation_gradient(3, 2)
    real(kind=kreal) :: current_jacobian, material_jacobian
 
    call shellmitc_directorcontributions(nn, director, zeta, shapefunc, shapederiv, director_contribution)
    call shellmitc_covariantbasis(nn, evaluation_coords, director, zeta, shapefunc, shapederiv, covariant_basis)
    if( abs(shellmitc_covariantjacobian(covariant_basis)) <= tiny(1.0d0) ) &
      stop "Invalid shell Jacobian"
    if( use_green_lagrange ) then
      translation_gradient = matmul(translation, shapederiv)
      covariant_basis(:, shell_xi:shell_eta) = covariant_basis(:, shell_xi:shell_eta)+translation_gradient
    endif
    call shellmitc_basisfromcovariant(covariant_basis, local_basis, reciprocal_basis, current_jacobian)
 
    tangent_basis = covariant_basis(:, shell_xi:shell_eta)
    material_local_basis = local_basis
    material_reciprocal_basis = reciprocal_basis
    integration_jacobian = current_jacobian
 
    if( .not. use_green_lagrange ) return
 
    call shellmitc_covariantbasis(nn, reference_coords, reference_director, zeta, shapefunc, shapederiv, reference_basis)
    call shellmitc_basisfromcovariant(reference_basis, material_local_basis, material_reciprocal_basis, material_jacobian)
    integration_jacobian = shellmitc_covariantjacobian(reference_basis)
 
    ! Finite-rotation MITC3/MITC4 use the current covariant tangent basis built above.
    if( etype == fe_mitc9_shell ) then
      tangent_basis = tangent_basis + translation_gradient
    else if( .not. fstr_is_finite_rotation_shell_element(etype, nn) ) then
      covariant_basis(:, shell_xi:shell_eta) = covariant_basis(:, shell_xi:shell_eta) + translation_gradient
      tangent_basis = covariant_basis(:, shell_xi:shell_eta)
    endif
  end subroutine shellmitc_preparepointkinematics
 
  !--------------------------------------------------------------------
  subroutine shellmitc_evaluatepointstrain(nn, zeta, shapefunc, shapederiv, translation, &
      director_increment, covariant_basis, use_green_lagrange, strain, reference_basis, &
      current_basis, reference_jacobian, current_jacobian)
    implicit none
 
    integer(kind=kint), intent(in) :: nn
    real(kind=kreal), intent(in) :: zeta, shapefunc(nn), shapederiv(nn, 2)
    real(kind=kreal), intent(in) :: translation(3, nn), director_increment(3, nn)
    real(kind=kreal), intent(in) :: covariant_basis(3, 3)
    logical, intent(in) :: use_green_lagrange
    real(kind=kreal), intent(out) :: strain(5)
    real(kind=kreal), intent(out) :: reference_basis(3, 3), current_basis(3, 3)
    real(kind=kreal), intent(out) :: reference_jacobian, current_jacobian
 
    real(kind=kreal) :: displacement_gradient(3, 3), translation_gradient(3, 2)
 
    displacement_gradient(:, shell_xi:shell_eta) = matmul(translation+zeta*director_increment, shapederiv)
    displacement_gradient(:, shell_zeta) = matmul(director_increment, shapefunc)
 
    reference_basis = covariant_basis
    current_basis = covariant_basis
 
    if( use_green_lagrange ) then
      translation_gradient = matmul(translation, shapederiv)
      current_basis(:, shell_xi:shell_eta) = covariant_basis(:, shell_xi:shell_eta)+translation_gradient
      reference_basis = current_basis-displacement_gradient
 
      strain = 0.5d0*shellmitc_symmetriccomponents( matmul(transpose(current_basis), current_basis) &
        -matmul(transpose(reference_basis), reference_basis))
    else
      strain = shellmitc_symmetriccomponents( matmul(transpose(covariant_basis), displacement_gradient))
    endif
 
    reference_jacobian = shellmitc_covariantjacobian(reference_basis)
    current_jacobian = shellmitc_covariantjacobian(current_basis)
 
  end subroutine shellmitc_evaluatepointstrain
 
  !--------------------------------------------------------------------
  !--------------------------------------------------------------------
  real(kind=kreal) function shellplanestresstracecoeff( gauss, n_layer )
    use mmechgauss
    use mmaterial, only: m_poisson, mc_isoelastic, getelastictype
    implicit none
 
    type(tgaussstatus), intent(in) :: gauss
    integer(kind=kint), intent(in) :: n_layer
 
    real(kind=kreal) :: nu, outa(2)
    logical :: ierr
 
    ! Isotropic plane-stress trace coefficient for the UL stress update.
    shellplanestresstracecoeff = 1.0d0
    if( .not. associated( gauss%pMaterial ) ) return
    if( getelastictype( gauss%pMaterial%mtype ) == 1 ) then
      stop "MITC shell UL orthotropic trace correction is not supported"
    endif
    nu = gauss%pMaterial%variables(m_poisson)
    call fetch_tabledata(mc_isoelastic, gauss%pMaterial%dict, outa, ierr)
    if( associated( gauss%pMaterial%shell_var ) ) then
      if( n_layer >= 1 .and. n_layer <= size( gauss%pMaterial%shell_var ) ) then
        if( gauss%pMaterial%shell_var(n_layer)%ortho == 0 ) then
          if( ierr ) then
            nu = gauss%pMaterial%shell_var(n_layer)%pp
          else
            nu = outa(2)
          endif
        else
          stop "MITC shell UL orthotropic trace correction is not supported"
        endif
      else if( .not. ierr ) then
        nu = outa(2)
      endif
    else if( .not. ierr ) then
      nu = outa(2)
    endif
 
    if( abs(1.0d0-nu) > 1.0d-12 ) then
      shellplanestresstracecoeff = (1.0d0-2.0d0*nu)/(1.0d0-nu)
    endif
 
  end function shellplanestresstracecoeff
 
 
  !--------------------------------------------------------------------
  subroutine shellmitc_setupreferencedirectors(etype, nn, thick, elem, v1, v2, v3, director)
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn
    real(kind=kreal), intent(in) :: thick, elem(3, nn)
    real(kind=kreal), intent(out) :: v1(3, nn), v2(3, nn), v3(3, nn), director(3, nn)
 
    integer :: nb
    real(kind=kreal) :: naturalcoord(2), nncoord(nn, 2), shapederiv(nn, 2)
    real(kind=kreal) :: g1(3), g2(3), e0(3), normv
 
    call getnodalnaturalcoord(etype, nncoord)
    naturalcoord = 0.0d0
    call getshapederiv(etype, naturalcoord, shapederiv)
    e0 = matmul(elem, shapederiv(:, 1))
 
    do nb = 1, nn
      naturalcoord = nncoord(nb, :)
      call getshapederiv(etype, naturalcoord, shapederiv)
      g1 = matmul(elem, shapederiv(:, 1))
      g2 = matmul(elem, shapederiv(:, 2))
 
      call cross_product(g1, g2, v3(:, nb))
      normv = dsqrt(dot_product(v3(:, nb), v3(:, nb)))
      v3(:, nb) = v3(:, nb)/normv
 
      call cross_product(v3(:, nb), e0, v2(:, nb))
      normv = dsqrt(dot_product(v2(:, nb), v2(:, nb)))
      if (normv > 1.0d-15) then
        v2(:, nb) = v2(:, nb)/normv
        call cross_product(v2(:, nb), v3(:, nb), v1(:, nb))
        normv = dsqrt(dot_product(v1(:, nb), v1(:, nb)))
        v1(:, nb) = v1(:, nb)/normv
      else
        v1(:, nb) = (/ 0.0d0, 0.0d0, -1.0d0 /)
        v2(:, nb) = (/ 0.0d0, 1.0d0, 0.0d0 /)
      endif
 
      call cross_product(v1(:, nb), v2(:, nb), v3(:, nb))
      normv = dsqrt(dot_product(v3(:, nb), v3(:, nb)))
      v3(:, nb) = v3(:, nb)/normv
      director(:, nb) = 0.5d0*thick*v3(:, nb)
    end do
  end subroutine shellmitc_setupreferencedirectors
 
  !--------------------------------------------------------------------
  subroutine shellmitc_buildinterpolationmatrix(nn, ndof, zeta, shapefunc, director, N)
    implicit none
 
    integer(kind=kint), intent(in) :: nn, ndof
    real(kind=kreal), intent(in) :: zeta, shapefunc(nn), director(3, nn)
    real(kind=kreal), intent(out) :: n(3, ndof*nn)
 
    integer :: nb, jsize
    real(kind=kreal) :: urot(3)
 
    n = 0.0d0
    do nb = 1, nn
      jsize = ndof*(nb-1)
      n(1, jsize+1) = shapefunc(nb)
      n(2, jsize+2) = shapefunc(nb)
      n(3, jsize+3) = shapefunc(nb)
      if (ndof >= 6) then
        urot = zeta*shapefunc(nb)*director(:, nb)
        ! ShellRotationGradient(v) = -ShellSkewMatrix(v) (transpose of a skew matrix).
        n(:, jsize+4:jsize+6) = -shellskewmatrix(urot)
      endif
    end do
  end subroutine shellmitc_buildinterpolationmatrix
 
 
  !--------------------------------------------------------------------
  subroutine shellmitc_setupnodaldirectors(etype, nn, thick, flag, elem, shell_disp, &
      finite_rotation_director, use_director_tangent, need_second_tangent, ndtriad, ndreftriad, ndcurtriad, &
      v1, v2, v3, a_over_2_v3, a_over_2_v3_ref, a_over_2_v3_deriv, a_over_2_v3_second)
    use mmaterial, only: updatelag
    implicit none
 
    integer, intent(in) :: flag
    integer(kind=kint), intent(in) :: etype, nn
    real(kind=kreal), intent(in) :: thick, elem(3, nn), shell_disp(6, nn)
    logical, intent(in) :: finite_rotation_director, use_director_tangent, need_second_tangent
    real(kind=kreal), intent(in), optional :: ndtriad(9, nn), ndreftriad(9, nn), ndcurtriad(9, nn)
    real(kind=kreal), intent(out) :: v1(3, nn), v2(3, nn), v3(3, nn)
    real(kind=kreal), intent(out) :: a_over_2_v3(3, nn), a_over_2_v3_ref(3, nn)
    real(kind=kreal), intent(out) :: a_over_2_v3_deriv(3, 3, nn), a_over_2_v3_second(3, 3, 3, nn)
 
    integer :: nb
    real(kind=kreal) :: director_ref(3), rotmat(3, 3)
    real(kind=kreal) :: nddirector(3, nn), ndrefdirector(3, nn), ndcurdirector(3, nn)
    logical :: has_nddirector, has_ndrefdirector, has_ndcurdirector
 
    has_nddirector = present(ndtriad)
    has_ndrefdirector = present(ndreftriad)
    has_ndcurdirector = present(ndcurtriad)
    if (has_nddirector)    nddirector(1:3, 1:nn)    = 0.5d0*thick*ndtriad(7:9, 1:nn)
    if (has_ndrefdirector) ndrefdirector(1:3, 1:nn) = 0.5d0*thick*ndreftriad(7:9, 1:nn)
    if (has_ndcurdirector) ndcurdirector(1:3, 1:nn) = 0.5d0*thick*ndcurtriad(7:9, 1:nn)
 
    a_over_2_v3_deriv = 0.0d0
    a_over_2_v3_second = 0.0d0
    call shellmitc_setupreferencedirectors(etype, nn, thick, elem, v1, v2, v3, a_over_2_v3)
 
    do nb = 1, nn
      a_over_2_v3_ref(:, nb) = a_over_2_v3(:, nb)
      if (finite_rotation_director .and. flag == updatelag .and. has_ndcurdirector) then
        a_over_2_v3_ref(:, nb) = ndcurdirector(:, nb)
      else if (finite_rotation_director .and. has_ndrefdirector) then
        a_over_2_v3_ref(:, nb) = ndrefdirector(:, nb)
      endif
 
      director_ref = a_over_2_v3_ref(:, nb)
      if (finite_rotation_director) then
        if (has_nddirector) then
          a_over_2_v3(:, nb) = nddirector(:, nb)
        else
          call shellrotationvectortomatrix(shell_disp(4:6, nb), rotmat)
          a_over_2_v3(:, nb) = matmul(rotmat, director_ref)
        endif
        if (use_director_tangent) then
          ! ShellRotationGradient(v) = -ShellSkewMatrix(v) (transpose of a skew matrix).
          a_over_2_v3_deriv(:, :, nb) = -shellskewmatrix(a_over_2_v3(:, nb))
          if( need_second_tangent ) then
            call shelldirectorincrementalsecondderiv(a_over_2_v3(:, nb), a_over_2_v3_second(:, :, :, nb))
          endif
        endif
      endif
 
      if (.not. use_director_tangent) then
        a_over_2_v3_deriv(:, :, nb) = -shellskewmatrix(a_over_2_v3(:, nb))
      endif
    end do
  end subroutine shellmitc_setupnodaldirectors
 
  !--------------------------------------------------------------------
  subroutine shellmitc_evaluatetyingpointvariation(etype, nn, ndof, tying_set, tying_point, &
      zeta_tying, elem, shell_disp, director, director_tangent, use_green_lagrange, &
      use_director_tangent, point_B, point_basis_variation, director_second_tangent, point_director_second_variation)
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn, ndof, tying_set, tying_point
    real(kind=kreal), intent(in) :: zeta_tying
    real(kind=kreal), intent(in) :: elem(3, nn), shell_disp(6, nn), director(3, nn)
    real(kind=kreal), intent(in) :: director_tangent(3, 3, nn)
    logical, intent(in) :: use_green_lagrange, use_director_tangent
    real(kind=kreal), intent(out) :: point_b(5, ndof*nn)
    real(kind=kreal), intent(out) :: point_basis_variation(3, ndof*nn, 3)
    real(kind=kreal), intent(in), optional :: director_second_tangent(3, 3, 3, nn)
    real(kind=kreal), intent(out), optional :: point_director_second_variation(5, 3, 3, nn)
 
    real(kind=kreal) :: naturalcoord(2), shapefunc(nn), shapederiv(nn, 2)
    real(kind=kreal) :: director_contribution(3, nn, 3)
    real(kind=kreal) :: kinematic_basis(3, 3), membrane_basis(3, 2)
    real(kind=kreal) :: translation_gradient(3, 2)
 
    call gettyingpoint(etype, tying_set, tying_point, naturalcoord)
    call getshapefunc(etype, naturalcoord, shapefunc)
    call getshapederiv(etype, naturalcoord, shapederiv)
 
    call shellmitc_directorcontributions(nn, director, zeta_tying, shapefunc, shapederiv, director_contribution)
    call shellmitc_covariantbasis(nn, elem, director, zeta_tying, shapefunc, shapederiv, kinematic_basis)
    if( abs(shellmitc_covariantjacobian(kinematic_basis)) <= tiny(1.0d0) ) &
      stop "Invalid shell Jacobian"
 
    membrane_basis = kinematic_basis(:, shell_xi:shell_eta)
    if( use_green_lagrange ) then
      translation_gradient = matmul(shell_disp(1:3, :), shapederiv)
      if( etype == fe_mitc9_shell ) then
        membrane_basis = membrane_basis+translation_gradient
      else
        kinematic_basis(:, shell_xi:shell_eta) = kinematic_basis(:, shell_xi:shell_eta)+translation_gradient
        membrane_basis = kinematic_basis(:, shell_xi:shell_eta)
      endif
    endif
 
    call shellmitc_buildfirststrainvariation(nn, ndof, zeta_tying, shapefunc, shapederiv, &
      kinematic_basis, membrane_basis, director_contribution, director_tangent, &
      use_director_tangent, point_b, point_basis_variation)
 
    if( use_director_tangent .and. present(director_second_tangent) .and. present(point_director_second_variation) ) then
      call shellmitc_builddirectorsecondvariation(nn, zeta_tying, shapefunc, shapederiv, &
        kinematic_basis, director_second_tangent, point_director_second_variation)
    endif
  end subroutine shellmitc_evaluatetyingpointvariation
 
  !--------------------------------------------------------------------
  subroutine shellmitc_evaluatetyingbatzeta(etype, nn, ndof, zeta_tying, elem, shell_disp, &
      director, director_tangent, use_green_lagrange, use_director_tangent, tying_B)
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn, ndof
    real(kind=kreal), intent(in) :: zeta_tying
    real(kind=kreal), intent(in) :: elem(3, nn), shell_disp(6, nn), director(3, nn)
    real(kind=kreal), intent(in) :: director_tangent(3, 3, nn)
    logical, intent(in) :: use_green_lagrange, use_director_tangent
    real(kind=kreal), intent(out) :: tying_b(5, ndof*nn, 6, 3)
 
    integer :: tying_set, tying_point
    real(kind=kreal) :: point_b(5, ndof*nn)
    real(kind=kreal) :: point_basis_variation(3, ndof*nn, 3)
 
    tying_b = 0.0d0
    do tying_set = 1, numoftyingsets(etype)
      do tying_point = 1, numoftyingpoints(etype, tying_set)
        call shellmitc_evaluatetyingpointvariation(etype, nn, ndof, tying_set, tying_point, &
          zeta_tying, elem, shell_disp, director, director_tangent, use_green_lagrange, &
          use_director_tangent, point_b, point_basis_variation)
        tying_b(:, :, tying_point, tying_set) = point_b
      end do
    end do
  end subroutine shellmitc_evaluatetyingbatzeta
 
  !--------------------------------------------------------------------
  subroutine shellmitc_evaluatetyingstiffnessdataatzeta(etype, nn, ndof, zeta_tying, &
      elem, shell_disp, director, director_tangent, director_second_tangent, &
      use_green_lagrange, use_director_tangent, tying_B, tying_basis_variation, tying_director_second_variation)
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn, ndof
    real(kind=kreal), intent(in) :: zeta_tying
    real(kind=kreal), intent(in) :: elem(3, nn), shell_disp(6, nn), director(3, nn)
    real(kind=kreal), intent(in) :: director_tangent(3, 3, nn)
    real(kind=kreal), intent(in) :: director_second_tangent(3, 3, 3, nn)
    logical, intent(in) :: use_green_lagrange, use_director_tangent
    real(kind=kreal), intent(out) :: tying_b(5, ndof*nn, 6, 3)
    real(kind=kreal), intent(out) :: tying_basis_variation(3, ndof*nn, 3, 6, 3)
    real(kind=kreal), intent(out) :: tying_director_second_variation(5, 3, 3, nn, 6, 3)
 
    integer :: tying_set, tying_point
    real(kind=kreal) :: point_b(5, ndof*nn)
    real(kind=kreal) :: point_basis_variation(3, ndof*nn, 3)
    real(kind=kreal) :: point_director_second_variation(5, 3, 3, nn)
 
    tying_b = 0.0d0
    tying_basis_variation = 0.0d0
    tying_director_second_variation = 0.0d0
 
    do tying_set = 1, numoftyingsets(etype)
      do tying_point = 1, numoftyingpoints(etype, tying_set)
        if( use_director_tangent ) then
          call shellmitc_evaluatetyingpointvariation(etype, nn, ndof, tying_set, tying_point, &
            zeta_tying, elem, shell_disp, director, director_tangent, use_green_lagrange, &
            use_director_tangent, point_b, point_basis_variation, director_second_tangent, point_director_second_variation)
          tying_director_second_variation(:, :, :, :, tying_point, tying_set) = point_director_second_variation
        else
          call shellmitc_evaluatetyingpointvariation(etype, nn, ndof, tying_set, tying_point, &
            zeta_tying, elem, shell_disp, director, director_tangent, use_green_lagrange, &
            use_director_tangent, point_b, point_basis_variation)
        endif
        tying_b(:, :, tying_point, tying_set) = point_b
        tying_basis_variation(:, :, :, tying_point, tying_set) = point_basis_variation
      end do
    end do
  end subroutine shellmitc_evaluatetyingstiffnessdataatzeta
 
  !--------------------------------------------------------------------
  subroutine shellmitc_buildfirststrainvariation(nn, ndof, zeta, shapefunc, shapederiv, &
      kinematic_basis, membrane_basis, director_contribution, director_tangent, &
      use_director_tangent, strain_variation, basis_variation)
    implicit none
 
    integer(kind=kint), intent(in) :: nn, ndof
    real(kind=kreal), intent(in) :: zeta, shapefunc(nn), shapederiv(nn, 2)
    real(kind=kreal), intent(in) :: kinematic_basis(3, 3), membrane_basis(3, 2)
    real(kind=kreal), intent(in) :: director_contribution(3, nn, 3)
    real(kind=kreal), intent(in) :: director_tangent(3, 3, nn)
    logical, intent(in) :: use_director_tangent
    real(kind=kreal), intent(out) :: strain_variation(5, ndof*nn)
    real(kind=kreal), intent(out) :: basis_variation(3, ndof*nn, 3)
 
    integer :: node, component, dof_index, rotation_start
 
    basis_variation = 0.0d0
    do node = 1, nn
      do component = 1, min(3, ndof)
        dof_index = ndof*(node-1)+component
        basis_variation(component, dof_index, shell_xi:shell_eta) = shapederiv(node, :)
      end do
      if( ndof < 6 ) cycle
 
      rotation_start = ndof*(node-1)+4
      if( use_director_tangent ) then
        basis_variation(:, rotation_start:rotation_start+2, shell_xi) = &
          shapederiv(node, shell_xi)*zeta*director_tangent(:, :, node)
        basis_variation(:, rotation_start:rotation_start+2, shell_eta) = &
          shapederiv(node, shell_eta)*zeta*director_tangent(:, :, node)
        basis_variation(:, rotation_start:rotation_start+2, shell_zeta) = shapefunc(node)*director_tangent(:, :, node)
      else
        ! ShellRotationGradient(v) = -ShellSkewMatrix(v) (transpose of a skew matrix).
        basis_variation(:, rotation_start:rotation_start+2, shell_xi) = &
          -shellskewmatrix(director_contribution(:, node, shell_xi))
        basis_variation(:, rotation_start:rotation_start+2, shell_eta) = &
          -shellskewmatrix(director_contribution(:, node, shell_eta))
        basis_variation(:, rotation_start:rotation_start+2, shell_zeta) = &
          -shellskewmatrix(director_contribution(:, node, shell_zeta))
      endif
    end do
 
    do dof_index = 1, ndof*nn
      strain_variation(:, dof_index) = shellmitc_symmetriccomponents( &
        matmul(transpose(basis_variation(:, dof_index, :)), kinematic_basis))
    end do
 
    ! In the TL MITC9 formulation only the translational membrane columns use
    ! the translated metric basis; rotational columns retain the director basis.
    do node = 1, nn
      do component = 1, min(3, ndof)
        dof_index = ndof*(node-1)+component
        strain_variation(1, dof_index) = shapederiv(node, shell_xi)*membrane_basis(component, shell_xi)
        strain_variation(2, dof_index) = shapederiv(node, shell_eta)*membrane_basis(component, shell_eta)
        strain_variation(3, dof_index) = shapederiv(node, shell_xi)*membrane_basis(component, shell_eta) &
          +shapederiv(node, shell_eta)*membrane_basis(component, shell_xi)
      end do
    end do
  end subroutine shellmitc_buildfirststrainvariation
 
  !--------------------------------------------------------------------
  subroutine shellmitc_builddirectorsecondvariation(nn, zeta, shapefunc, shapederiv, &
      kinematic_basis, director_second_tangent, director_second_variation)
    implicit none
 
    integer(kind=kint), intent(in) :: nn
    real(kind=kreal), intent(in) :: zeta, shapefunc(nn), shapederiv(nn, 2)
    real(kind=kreal), intent(in) :: kinematic_basis(3, 3)
    real(kind=kreal), intent(in) :: director_second_tangent(3, 3, 3, nn)
    real(kind=kreal), intent(out) :: director_second_variation(5, 3, 3, nn)
 
    integer :: node, m, n
    real(kind=kreal) :: second_basis(3, 3)
 
    director_second_variation = 0.0d0
    do node = 1, nn
      do n = 1, 3
        do m = 1, 3
          second_basis(:, shell_xi) = zeta*shapederiv(node, shell_xi) *director_second_tangent(:, m, n, node)
          second_basis(:, shell_eta) = zeta*shapederiv(node, shell_eta) *director_second_tangent(:, m, n, node)
          second_basis(:, shell_zeta) = shapefunc(node) *director_second_tangent(:, m, n, node)
          director_second_variation(:, m, n, node) = shellmitc_symmetriccomponents( &
            matmul(transpose(second_basis), kinematic_basis))
        end do
      end do
    end do
  end subroutine shellmitc_builddirectorsecondvariation
 
  !--------------------------------------------------------------------
  subroutine shellmitc_evaluatetyingoperator(etype, xi, eta, nterms, target_component, &
      source_component, tying_point, tying_group, coefficient, replaced_component)
    implicit none
 
    integer(kind=kint), intent(in) :: etype
    real(kind=kreal), intent(in) :: xi, eta
    integer, intent(out) :: nterms
    integer, intent(out) :: target_component(max_tying_terms), source_component(MAX_TYING_TERMS)
    integer, intent(out) :: tying_point(MAX_TYING_TERMS), tying_group(MAX_TYING_TERMS)
    real(kind=kreal), intent(out) :: coefficient(max_tying_terms)
    logical, intent(out) :: replaced_component(5)
 
    integer :: ip, k
    real(kind=kreal) :: xxi, eeta, h
 
    nterms = 0
 
    if( etype == fe_mitc4_shell ) then
      nterms = 4
      target_component(1:4) = (/ 4, 4, 5, 5 /)
      source_component(1:4) = (/ 4, 4, 5, 5 /)
      tying_point(1:4) = (/ 4, 2, 1, 3 /)
      tying_group(1:4) = (/ 1, 1, 1, 1 /)
      coefficient(1:4) = (/ 0.5d0*( 1.0d0-xi ), 0.5d0*( 1.0d0+xi ), 0.5d0*( 1.0d0-eta ), 0.5d0*( 1.0d0+eta ) /)
 
    else if( etype == fe_mitc9_shell ) then
      xxi  = xi /dsqrt( 1.0d0/3.0d0 )
      eeta = eta/dsqrt( 3.0d0/5.0d0 )
      do ip = 1, numoftyingpoints(etype, 1)
        h = ( 0.5d0*( 1.0d0+mitc9_xi_sign(ip, 1)*xxi ) )   &
          *( ( 0.5d0*mitc9_eta_sign(ip, 1)*eeta ) *( 1.0d0+mitc9_eta_sign(ip, 1)*eeta )       &
          +( 1.0d0-mitc9_eta_sign(ip, 1)*mitc9_eta_sign(ip, 1) ) *( 1.0d0-eeta*eeta ) )
        target_component(nterms+1:nterms+2) = (/ 1, 5 /)
        source_component(nterms+1:nterms+2) = (/ 1, 5 /)
        tying_point(nterms+1:nterms+2) = ip
        tying_group(nterms+1:nterms+2) = 1
        coefficient(nterms+1:nterms+2) = h
        nterms = nterms+2
      end do
 
      xxi  = xi /dsqrt( 3.0d0/5.0d0 )
      eeta = eta/dsqrt( 1.0d0/3.0d0 )
      do ip = 1, numoftyingpoints(etype, 2)
        h = ( ( 0.5d0*mitc9_xi_sign(ip, 2)*xxi ) *( 1.0d0+mitc9_xi_sign(ip, 2)*xxi ) &
          +( 1.0d0-mitc9_xi_sign(ip, 2)*mitc9_xi_sign(ip, 2) ) *( 1.0d0-xxi*xxi ) ) &
          *( 0.5d0*( 1.0d0+mitc9_eta_sign(ip, 2)*eeta ) )
        target_component(nterms+1:nterms+2) = (/ 2, 4 /)
        source_component(nterms+1:nterms+2) = (/ 2, 4 /)
        tying_point(nterms+1:nterms+2) = ip
        tying_group(nterms+1:nterms+2) = 2
        coefficient(nterms+1:nterms+2) = h
        nterms = nterms+2
      end do
 
      xxi  = xi /dsqrt( 1.0d0/3.0d0 )
      eeta = eta/dsqrt( 1.0d0/3.0d0 )
      do ip = 1, numoftyingpoints(etype, 3)
        h = ( 0.5d0*( 1.0d0+mitc9_xi_sign(ip, 1)*xxi ) ) *( 0.5d0*( 1.0d0+mitc9_eta_sign(ip, 1)*eeta ) )
        nterms = nterms+1
        target_component(nterms) = 3
        source_component(nterms) = 3
        tying_point(nterms) = ip
        tying_group(nterms) = 3
        coefficient(nterms) = h
      end do
 
    else if( etype == fe_mitc3_shell ) then
      nterms = 8
      target_component(1:8) = (/ 4, 4, 4, 4, 5, 5, 5, 5 /)
      source_component(1:8) = (/ 4, 5, 4, 5, 4, 5, 4, 5 /)
      tying_point(1:8) = (/ 2, 1, 3, 3, 2, 1, 3, 3 /)
      tying_group(1:8) = (/ 1, 1, 1, 1, 1, 1, 1, 1 /)
      coefficient(1:8) = (/ 1.0d0-xi, xi, xi, -xi, eta, 1.0d0-eta, -eta, eta /)
    endif
 
    replaced_component(:) = .false.
    do k = 1, nterms
      replaced_component(target_component(k)) = .true.
    end do
 
  end subroutine shellmitc_evaluatetyingoperator
 
  !--------------------------------------------------------------------
  subroutine shellmitc_applyassumedstrain(etype, nn, ndof, xi_lx, eta_lx, tying_B, B, &
      tying_director_second_variation, B2rot)
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn, ndof
    real(kind=kreal), intent(in) :: xi_lx, eta_lx
    real(kind=kreal), intent(in) :: tying_b(5, ndof*nn, 6, 3)
    real(kind=kreal), intent(inout) :: b(5, ndof*nn)
    real(kind=kreal), intent(in), optional :: tying_director_second_variation(5, 3, 3, nn, 6, 3)
    real(kind=kreal), intent(inout), optional :: b2rot(:, :, :, :)
 
    integer :: k, c, jsize, m, n, nb, nterms
    integer :: target_component(MAX_TYING_TERMS), source_component(MAX_TYING_TERMS)
    integer :: tying_point(MAX_TYING_TERMS), tying_group(MAX_TYING_TERMS)
    real(kind=kreal) :: coefficient(max_tying_terms)
    logical :: replaced_component(5)
 
    call shellmitc_evaluatetyingoperator(etype, xi_lx, eta_lx, nterms, target_component, &
      source_component, tying_point, tying_group, coefficient, replaced_component)
    if( nterms == 0 ) return
 
    do c = 1, 5
      if( replaced_component(c) ) b(c, 1:ndof*nn) = 0.0d0
    end do
    do k = 1, nterms
      do jsize = 1, ndof*nn
        b(target_component(k), jsize) = b(target_component(k), jsize) &
          +coefficient(k)*tying_b(source_component(k), jsize, tying_point(k), tying_group(k))
      end do
    end do
 
    if( present(tying_director_second_variation) .and. present(b2rot) ) then
      do c = 1, 5
        if( replaced_component(c) ) b2rot(c, :, :, :) = 0.0d0
      end do
      do k = 1, nterms
        do nb = 1, nn
          do n = 1, 3
            do m = 1, 3
              b2rot(target_component(k), m, n, nb) = b2rot(target_component(k), m, n, nb) &
                +coefficient(k)*tying_director_second_variation(source_component(k), m, n, nb, &
                tying_point(k), tying_group(k))
            end do
          end do
        end do
      end do
    endif
 
  end subroutine shellmitc_applyassumedstrain
 
  !--------------------------------------------------------------------
  pure real(kind=kreal) function shellmitc_strainsecondvariationpair(component, basis_variation_i, basis_variation_j)
    implicit none
 
    integer, intent(in) :: component
    real(kind=kreal), intent(in) :: basis_variation_i(3, 3), basis_variation_j(3, 3)
 
    select case( component )
      case( 1 )
        shellmitc_strainsecondvariationpair = dot_product( basis_variation_i(:, shell_xi), basis_variation_j(:, shell_xi))
      case( 2 )
        shellmitc_strainsecondvariationpair = dot_product( basis_variation_i(:, shell_eta), basis_variation_j(:, shell_eta))
      case( 3 )
        shellmitc_strainsecondvariationpair = dot_product( &
          basis_variation_i(:, shell_xi), basis_variation_j(:, shell_eta)) &
          +dot_product(basis_variation_i(:, shell_eta), basis_variation_j(:, shell_xi))
      case( 4 )
        shellmitc_strainsecondvariationpair = dot_product( &
          basis_variation_i(:, shell_eta), basis_variation_j(:, shell_zeta)) &
          +dot_product(basis_variation_i(:, shell_zeta), basis_variation_j(:, shell_eta))
      case( 5 )
        shellmitc_strainsecondvariationpair = dot_product( &
          basis_variation_i(:, shell_zeta), basis_variation_j(:, shell_xi)) &
          +dot_product(basis_variation_i(:, shell_xi), basis_variation_j(:, shell_zeta))
      case default
        shellmitc_strainsecondvariationpair = 0.0d0
    end select
 
  end function shellmitc_strainsecondvariationpair
 
  !--------------------------------------------------------------------
  subroutine shellmitc_addgeometricstiffness(etype, nn, ndof, kinematics, &
      use_director_tangent, use_green_lagrange, add_geometric_stiffness, &
      xi, eta, integration_weight, layer_weight, stress, B, basis_variation, &
      tying_basis_variation, reciprocal_basis, material_local_basis, material_reciprocal_basis, B2rot, director, stiff)
    use mmaterial, only: updatelag
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn, ndof, kinematics
    logical, intent(in) :: use_director_tangent, use_green_lagrange
    logical, intent(in) :: add_geometric_stiffness
    real(kind=kreal), intent(in) :: xi, eta, integration_weight, layer_weight
    real(kind=kreal), intent(in) :: stress(5), b(5, ndof*nn)
    real(kind=kreal), intent(in) :: basis_variation(3, ndof*nn, 3)
    real(kind=kreal), intent(in) :: tying_basis_variation(3, ndof*nn, 3, 6, 3)
    real(kind=kreal), intent(in) :: reciprocal_basis(3, 3)
    real(kind=kreal), intent(in) :: material_local_basis(3, 3)
    real(kind=kreal), intent(in) :: material_reciprocal_basis(3, 3)
    real(kind=kreal), intent(in) :: b2rot(5, 3, 3, nn)
    real(kind=kreal), intent(in) :: director(3, nn)
    real(kind=kreal), intent(inout) :: stiff(ndof*nn, ndof*nn)
 
    integer :: isize, jsize, k, c, ip, it, nterms
    integer :: target_component(MAX_TYING_TERMS), source_component(MAX_TYING_TERMS)
    integer :: tying_point(MAX_TYING_TERMS), tying_group(MAX_TYING_TERMS)
    real(kind=kreal) :: coefficient(max_tying_terms)
    logical :: replaced_component(5)
    real(kind=kreal) :: qf_stress_integrand(ndof*nn), vol_deriv(ndof*nn)
    real(kind=kreal) :: geo_term, local_stress(3, 3), s_global(3, 3)
    real(kind=kreal) :: basis_gradient(3, 3, ndof*nn)
    real(kind=kreal) :: stress_basis_gradient(3, 3, ndof*nn)
 
    if( kinematics == updatelag ) then
      do jsize = 1, ndof*nn
        vol_deriv(jsize) = sum(reciprocal_basis*basis_variation(:, jsize, :))
        qf_stress_integrand(jsize) = dot_product(stress, b(:, jsize))
      end do
    endif
 
    if( add_geometric_stiffness ) then
      if( use_green_lagrange .or. kinematics == updatelag ) then
        call shellmitc_evaluatetyingoperator(etype, xi, eta, nterms, target_component, &
          source_component, tying_point, tying_group, coefficient, replaced_component)
        do jsize = 1, ndof*nn
          do isize = 1, ndof*nn
            geo_term = 0.0d0
            do c = 1, 5
              if( replaced_component(c) ) cycle
              geo_term = geo_term+stress(c)*shellmitc_strainsecondvariationpair(c, &
                basis_variation(:, isize, :), basis_variation(:, jsize, :))
            end do
            do k = 1, nterms
              ip = tying_point(k)
              it = tying_group(k)
              geo_term = geo_term+stress(target_component(k))*coefficient(k) &
                *shellmitc_strainsecondvariationpair(source_component(k), tying_basis_variation(:, isize, :, ip, it), &
                tying_basis_variation(:, jsize, :, ip, it))
            end do
            stiff(isize, jsize) = stiff(isize, jsize) +integration_weight*layer_weight*geo_term
          end do
        end do
      else
        local_stress = reshape((/ stress(1), stress(3), stress(5), &
          stress(3), stress(2), stress(4), stress(5), stress(4), 0.0d0 /), (/ 3, 3 /))
        s_global = matmul(material_local_basis, matmul(local_stress, transpose(material_local_basis)))
 
        do jsize = 1, ndof*nn
          basis_gradient(:, :, jsize) = matmul(basis_variation(:, jsize, :), transpose(material_reciprocal_basis))
          stress_basis_gradient(:, :, jsize) = matmul(basis_gradient(:, :, jsize), s_global)
        end do
        do jsize = 1, ndof*nn
          do isize = 1, ndof*nn
            stiff(isize, jsize) = stiff(isize, jsize) +integration_weight*layer_weight*sum( &
              basis_gradient(:, :, isize)*stress_basis_gradient(:, :, jsize))
          end do
        end do
      endif
      if( kinematics == updatelag ) then
        do jsize = 1, ndof*nn
          do isize = 1, ndof*nn
            stiff(isize, jsize) = stiff(isize, jsize) &
              +integration_weight*layer_weight*qf_stress_integrand(isize)*vol_deriv(jsize)
          end do
        end do
      endif
      if( use_director_tangent .and. (use_green_lagrange .or. kinematics == updatelag) ) then
        call shellmitc_adddirectorstressstiffness(nn, ndof, integration_weight, layer_weight, &
          stress, b2rot, director, stiff)
      endif
    endif
 
  end subroutine shellmitc_addgeometricstiffness
 
  !--------------------------------------------------------------------
  subroutine shellmitc_builddrillingsecondvariation(nn, zeta, shapefunc, shapederiv, &
      director_second_tangent, reciprocal_basis, local_basis, use_director_tangent, drilling_second_variation)
    implicit none
 
    integer(kind=kint), intent(in) :: nn
    real(kind=kreal), intent(in) :: zeta, shapefunc(nn), shapederiv(nn, 2)
    real(kind=kreal), intent(in) :: director_second_tangent(3, 3, 3, nn)
    real(kind=kreal), intent(in) :: reciprocal_basis(3, 3), local_basis(3, 3)
    logical, intent(in) :: use_director_tangent
    real(kind=kreal), intent(out) :: drilling_second_variation(3, 3, nn)
 
    integer :: node, m, n
    real(kind=kreal) :: second_basis(3, 3), skew_gradient(3, 3)
 
    drilling_second_variation = 0.0d0
    if( .not. use_director_tangent ) return
    do node = 1, nn
      do n = 1, 3
        do m = 1, 3
          second_basis(:, shell_xi) = shapederiv(node, shell_xi)*zeta*director_second_tangent(:, m, n, node)
          second_basis(:, shell_eta) = shapederiv(node, shell_eta)*zeta*director_second_tangent(:, m, n, node)
          second_basis(:, shell_zeta) = shapefunc(node)*director_second_tangent(:, m, n, node)
          skew_gradient = matmul(reciprocal_basis, transpose(second_basis))
          skew_gradient = skew_gradient-transpose(skew_gradient)
          drilling_second_variation(m, n, node) = dot_product( local_basis(:, shell_xi), &
            matmul(skew_gradient, local_basis(:, shell_eta)))
        end do
      end do
    end do
 
  end subroutine shellmitc_builddrillingsecondvariation
 
  !--------------------------------------------------------------------
  subroutine shellmitc_builddrillingvector(nn, ndof, finite_rotation_director, shapefunc, &
      point_triad, reciprocal_basis, basis_variation, displacement, director, nddrill, Cv, Cv_disp)
    implicit none
 
    integer(kind=kint), intent(in) :: nn, ndof
    logical, intent(in) :: finite_rotation_director
    real(kind=kreal), intent(in) :: shapefunc(nn), point_triad(3, 3)
    real(kind=kreal), intent(in) :: reciprocal_basis(3, 3)
    real(kind=kreal), intent(in) :: basis_variation(3, ndof*nn, 3)
    real(kind=kreal), intent(in) :: displacement(ndof*nn), director(3, nn)
    real(kind=kreal), intent(in), optional :: nddrill(nn)
    real(kind=kreal), intent(out) :: cv(ndof*nn), cv_disp
 
    integer :: j, nb, jrot
    real(kind=kreal) :: cv_w(ndof*nn), cv_theta(ndof*nn)
    real(kind=kreal) :: cmat(3, 3), drill_axis(3), axis_norm
 
    do j = 1, ndof*nn
      cmat = matmul(reciprocal_basis, transpose(basis_variation(:, j, :)))
      cmat = cmat-transpose(cmat)
      cv_w(j) = dot_product(point_triad(:, shell_xi), matmul(cmat, point_triad(:, shell_eta)))
    end do
 
    cv_theta = 0.0d0
    if( ndof >= 6 ) then
      do nb = 1, nn
        jrot = ndof*(nb-1)+4
        drill_axis = point_triad(:, shell_zeta)
        if( finite_rotation_director .and. present(nddrill) ) then
          drill_axis = director(:, nb)
          axis_norm = sqrt(dot_product(drill_axis, drill_axis))
          if( axis_norm > 0.0d0 ) then
            drill_axis = drill_axis/axis_norm
          else
            drill_axis = point_triad(:, shell_zeta)
          endif
        endif
        cv_theta(jrot:jrot+2) = shapefunc(nb)*drill_axis
      end do
    endif
 
    cv = cv_theta-0.5d0*cv_w
    cv_disp = dot_product(cv, displacement)
    if( finite_rotation_director .and. present(nddrill) ) then
      cv_disp = -0.5d0*dot_product(cv_w, displacement)+dot_product(shapefunc, nddrill)
    endif
  end subroutine shellmitc_builddrillingvector
 
  !--------------------------------------------------------------------
  subroutine shellmitc_adddrillingstiffness(nn, ndof, finite_rotation_director, &
      use_director_tangent, shapefunc, Cv, Cv_disp, Cv_w_second, displacement, &
      director, director_deriv, alpha, integration_weight, layer_weight, nddrill, stiff)
    implicit none
 
    integer(kind=kint), intent(in) :: nn, ndof
    logical, intent(in) :: finite_rotation_director, use_director_tangent
    real(kind=kreal), intent(in) :: shapefunc(nn), cv(ndof*nn), cv_disp
    real(kind=kreal), intent(in) :: cv_w_second(3, 3, nn), displacement(ndof*nn)
    real(kind=kreal), intent(in) :: director(3, nn), director_deriv(3, 3, nn)
    real(kind=kreal), intent(in) :: alpha, integration_weight, layer_weight
    real(kind=kreal), intent(in), optional :: nddrill(nn)
    real(kind=kreal), intent(inout) :: stiff(ndof*nn, ndof*nn)
 
    integer :: isize, jsize, nb, m, n
    real(kind=kreal) :: scale, cv_deriv, cv_deriv_disp
    real(kind=kreal) :: drill_axis(3), drill_coeff, axis_norm
 
    scale = integration_weight*layer_weight*alpha
    do jsize = 1, ndof*nn
      do isize = 1, ndof*nn
        stiff(isize, jsize) = stiff(isize, jsize)+scale*cv(isize)*cv(jsize)
      end do
    end do
    if( .not. use_director_tangent ) return
 
    do nb = 1, nn
      do n = 1, 3
        jsize = ndof*(nb-1)+3+n
        cv_deriv_disp = 0.0d0
        do m = 1, 3
          isize = ndof*(nb-1)+3+m
          cv_deriv = -0.5d0*cv_w_second(m, n, nb)
          ! Cv_disp takes the drilling angle from nddrill, so the drill-axis rotation below enters only Cv.
          cv_deriv_disp = cv_deriv_disp+cv_deriv*displacement(isize)
          if( finite_rotation_director .and. present(nddrill) ) then
            drill_axis = director(:, nb)
            axis_norm = sqrt(dot_product(drill_axis, drill_axis))
            if( axis_norm > 0.0d0 ) then
              drill_axis = drill_axis/axis_norm
              drill_coeff = dot_product(drill_axis, director_deriv(:, n, nb))
              cv_deriv = cv_deriv+shapefunc(nb) *(director_deriv(m, n, nb)-drill_axis(m)*drill_coeff)/axis_norm
            endif
          endif
          stiff(isize, jsize) = stiff(isize, jsize)+scale*cv_deriv*cv_disp
        end do
        stiff(:, jsize) = stiff(:, jsize)+scale*cv*cv_deriv_disp
      end do
    end do
  end subroutine shellmitc_adddrillingstiffness
 
  !--------------------------------------------------------------------
  subroutine shellmitc_adddirectorstressstiffness(nn, ndof, w_w_w_det, layer_weight, Sv_force, B2rot, a_over_2_v3, stiff)
    implicit none
 
    integer(kind=kint), intent(in) :: nn, ndof
    real(kind=kreal), intent(in) :: w_w_w_det, layer_weight, sv_force(5)
    real(kind=kreal), intent(in) :: b2rot(5, 3, 3, nn), a_over_2_v3(3, nn)
    real(kind=kreal), intent(inout) :: stiff(ndof*nn, ndof*nn)
 
    integer :: i, j, m, n, nb, isize, jsize
    real(kind=kreal) :: drill_axis(3), rot_projector(3, 3), hrot(3, 3)
    real(kind=kreal) :: hess_coeff, drill_coeff, axis_norm
 
    do nb = 1, nn
      drill_axis(1:3) = a_over_2_v3(1:3, nb)
      axis_norm = sqrt(sum(drill_axis(1:3)*drill_axis(1:3)))
      if( axis_norm > 0.0d0 ) then
        drill_axis(1:3) = drill_axis(1:3)/axis_norm
      else
        drill_axis(1:3) = (/ 0.0d0, 0.0d0, 1.0d0 /)
      endif
      rot_projector(:, :) = 0.0d0
      do i = 1, 3
        rot_projector(i, i) = 1.0d0
      end do
      do i = 1, 3
        do j = 1, 3
          rot_projector(i, j) = rot_projector(i, j)-drill_axis(i)*drill_axis(j)
        end do
      end do
      do n = 1, 3
        do m = 1, 3
          hrot(m, n) = dot_product(sv_force(1:5), b2rot(1:5, m, n, nb))
        end do
      end do
      do n = 1, 3
        jsize = ndof*(nb-1)+3+n
        do m = 1, 3
          isize = ndof*(nb-1)+3+m
          hess_coeff = 0.0d0
          do j = 1, 3
            do i = 1, 3
              hess_coeff = hess_coeff +rot_projector(i, m)*hrot(i, j)*rot_projector(j, n)
            end do
          end do
          drill_coeff = 0.0d0
          do i = 1, 3
            drill_coeff = drill_coeff+drill_axis(i)*hrot(i, n)
          end do
          do j = 1, 3
            do i = 1, 3
              drill_coeff = drill_coeff +rot_projector(i, n)*hrot(i, j)*drill_axis(j)
            end do
          end do
          stiff(isize, jsize) = stiff(isize, jsize) +w_w_w_det*layer_weight *(hess_coeff+drill_axis(m)*drill_coeff)
        end do
      end do
    end do
 
  end subroutine shellmitc_adddirectorstressstiffness
 
  !--------------------------------------------------------------------
  subroutine shellmixeddofmap(mixflag, nn, ndof, sstable, is_mixed)
    implicit none
 
    integer(kind=kint), intent(in) :: mixflag, nn, ndof
    integer, intent(out) :: sstable(ndof*nn)
    logical, intent(out) :: is_mixed
 
    integer :: i
 
    is_mixed = .true.
    if( mixflag == 1 .and. ndof*nn == 24 ) then
      sstable = (/ 1, 2, 3, 7, 8, 9, 13, 14, 15, 19, 20, 21, 4, 5, 6, 10, 11, 12, 16, 17, 18, 22, 23, 24 /)
    else if( mixflag == 2 .and. ndof*nn == 18 ) then
      sstable = (/ 1, 2, 3, 7, 8, 9, 13, 14, 15, 4, 5, 6, 10, 11, 12, 16, 17, 18 /)
    else
      is_mixed = .false.
      sstable = [(i, i=1, ndof*nn)]
    endif
  end subroutine shellmixeddofmap
 
  !--------------------------------------------------------------------
  subroutine shellapplymixeddofordering(mixflag, nn, ndof, stiff, qf)
    implicit none
 
    integer(kind=kint), intent(in) :: mixflag, nn, ndof
    real(kind=kreal), intent(inout), optional :: stiff(ndof*nn, ndof*nn), qf(ndof*nn)
 
    integer :: i, j, sstable(ndof*nn)
    real(kind=kreal) :: stiff_work(ndof*nn, ndof*nn), qf_work(ndof*nn)
    logical :: is_mixed
 
    call shellmixeddofmap(mixflag, nn, ndof, sstable, is_mixed)
    if( .not. is_mixed ) return
 
    if( present(stiff) ) then
      stiff_work = stiff
      do j = 1, ndof*nn
        do i = 1, ndof*nn
          stiff(i, j) = stiff_work(sstable(i), sstable(j))
        end do
      end do
    endif
    if( present(qf) ) then
      qf_work = qf
      do i = 1, ndof*nn
        qf(i) = qf_work(sstable(i))
      end do
    endif
  end subroutine shellapplymixeddofordering
 
  !--------------------------------------------------------------------
  subroutine shellmitc_evaluatetyingstrains(etype, nn, zeta, elem, edisp, director, &
      director_increment, use_green_lagrange, tying_strain)
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn
    real(kind=kreal), intent(in) :: zeta
    real(kind=kreal), intent(in) :: elem(3, nn), edisp(6, nn)
    real(kind=kreal), intent(in) :: director(3, nn), director_increment(3, nn)
    logical, intent(in) :: use_green_lagrange
    real(kind=kreal), intent(out) :: tying_strain(5, 6, 3)
 
    integer :: tying_set, tying_point
    real(kind=kreal) :: tying_zeta, naturalcoord(2)
    real(kind=kreal) :: shapefunc(nn), shapederiv(nn, 2)
    real(kind=kreal) :: covariant_basis(3, 3), reference_basis(3, 3), current_basis(3, 3)
    real(kind=kreal) :: reference_jacobian, current_jacobian, point_strain(5)
 
    tying_strain = 0.0d0
    tying_zeta = shellmitc_tyingzeta(etype, zeta)
 
    do tying_set = 1, numoftyingsets(etype)
      do tying_point = 1, numoftyingpoints(etype, tying_set)
        call gettyingpoint(etype, tying_set, tying_point, naturalcoord)
        call getshapefunc(etype, naturalcoord, shapefunc)
        call getshapederiv(etype, naturalcoord, shapederiv)
        call shellmitc_covariantbasis(nn, elem, director, tying_zeta, shapefunc, shapederiv, covariant_basis)
        if( abs(shellmitc_covariantjacobian(covariant_basis)) <= tiny(1.0d0) ) &
          stop "Invalid shell Jacobian"
        call shellmitc_evaluatepointstrain(nn, tying_zeta, shapefunc, shapederiv, &
          edisp(1:3, :), director_increment, covariant_basis, use_green_lagrange, &
          point_strain, reference_basis, current_basis, reference_jacobian, current_jacobian)
        tying_strain(:, tying_point, tying_set) = point_strain
      end do
    end do
 
  end subroutine shellmitc_evaluatetyingstrains
 
  !--------------------------------------------------------------------
  subroutine shellmitc_setupstresskinematics(nn, ecoord, kinematics, finite_rotation, edisp, &
      nddirector, ndrefdirector, nddirector_deriv, ndbase_disp, elem, director, director_increment)
    use mmaterial, only: updatelag
    implicit none
 
    integer(kind=kint), intent(in) :: nn, kinematics
    logical, intent(in) :: finite_rotation
    real(kind=kreal), intent(in) :: ecoord(3, nn), edisp(6, nn)
    real(kind=kreal), intent(in) :: nddirector(3, nn), ndrefdirector(3, nn)
    real(kind=kreal), intent(in) :: nddirector_deriv(3, 3, nn)
    real(kind=kreal), intent(in), optional :: ndbase_disp(6, nn)
    real(kind=kreal), intent(out) :: elem(3, nn)
    real(kind=kreal), intent(out) :: director(3, nn), director_increment(3, nn)
 
    integer :: na
 
    elem = ecoord
    if( kinematics == updatelag .and. present(ndbase_disp) ) then
      elem = elem+ndbase_disp(1:3, :)
    endif
    if( kinematics == updatelag ) elem = elem+0.5d0*edisp(1:3, :)
 
    do na = 1, nn
      director(:, na) = nddirector(:, na)
      director_increment(:, na) = director(:, na)-ndrefdirector(:, na)
      if( .not. finite_rotation ) then
        director_increment(:, na) = matmul(nddirector_deriv(:, :, na), edisp(4:6, na))
      else if( kinematics == updatelag ) then
        director(:, na) = ndrefdirector(:, na)+0.5d0*director_increment(:, na)
      endif
    end do
  end subroutine shellmitc_setupstresskinematics
 
  !--------------------------------------------------------------------
  subroutine shellmitc_evaluateassumedstrain(etype, nn, naturalcoord, zeta, elem, edisp, &
      director, director_increment, use_green_lagrange, tying_strain, dstrain, &
      strain_tensor, covariant_basis, reciprocal_basis, material_local_basis, &
      material_reciprocal_basis, reference_basis, current_basis, reference_jacobian, current_jacobian)
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn
    real(kind=kreal), intent(in) :: naturalcoord(2), zeta
    real(kind=kreal), intent(in) :: elem(3, nn), edisp(6, nn)
    real(kind=kreal), intent(in) :: director(3, nn), director_increment(3, nn)
    logical, intent(in) :: use_green_lagrange
    real(kind=kreal), intent(in) :: tying_strain(5, 6, 3)
    real(kind=kreal), intent(out) :: dstrain(6), strain_tensor(3, 3)
    real(kind=kreal), intent(out) :: covariant_basis(3, 3), reciprocal_basis(3, 3)
    real(kind=kreal), intent(out) :: material_local_basis(3, 3)
    real(kind=kreal), intent(out) :: material_reciprocal_basis(3, 3)
    real(kind=kreal), intent(out) :: reference_basis(3, 3), current_basis(3, 3)
    real(kind=kreal), intent(out) :: reference_jacobian, current_jacobian
 
    integer :: k, c, nterms
    integer :: target_component(MAX_TYING_TERMS), source_component(MAX_TYING_TERMS)
    integer :: tying_point(MAX_TYING_TERMS), tying_group(MAX_TYING_TERMS)
    real(kind=kreal) :: coefficient(max_tying_terms)
    logical :: replaced_component(5)
    real(kind=kreal) :: shapefunc(nn), shapederiv(nn, 2), strain(5)
    real(kind=kreal) :: local_basis(3, 3), jacobian, material_jacobian
 
    call getshapefunc(etype, naturalcoord, shapefunc)
    call getshapederiv(etype, naturalcoord, shapederiv)
    call shellmitc_covariantbasis(nn, elem, director, zeta, shapefunc, shapederiv, covariant_basis)
    if( abs(shellmitc_covariantjacobian(covariant_basis)) <= tiny(1.0d0) ) &
      stop "Invalid shell Jacobian"
    call shellmitc_basisfromcovariant(covariant_basis, local_basis, reciprocal_basis, jacobian)
    call shellmitc_evaluatepointstrain(nn, zeta, shapefunc, shapederiv, edisp(1:3, :), &
      director_increment, covariant_basis, use_green_lagrange, strain, reference_basis, &
      current_basis, reference_jacobian, current_jacobian)
    call shellmitc_evaluatetyingoperator(etype, naturalcoord(1), naturalcoord(2), nterms, &
      target_component, source_component, tying_point, tying_group, coefficient, replaced_component)
    do c = 1, 5
      if( replaced_component(c) ) strain(c) = 0.0d0
    end do
    do k = 1, nterms
      strain(target_component(k)) = strain(target_component(k)) &
        +coefficient(k)*tying_strain(source_component(k), tying_point(k), tying_group(k))
    end do
 
    strain_tensor = 0.0d0
    strain_tensor(1, 1) = strain(1)
    strain_tensor(2, 2) = strain(2)
    strain_tensor(1, 2) = 0.5d0*strain(3)
    strain_tensor(2, 1) = strain_tensor(1, 2)
    strain_tensor(2, 3) = 0.5d0*strain(4)
    strain_tensor(3, 2) = strain_tensor(2, 3)
    strain_tensor(3, 1) = 0.5d0*strain(5)
    strain_tensor(1, 3) = strain_tensor(3, 1)
 
    material_local_basis = local_basis
    material_reciprocal_basis = reciprocal_basis
    if( use_green_lagrange ) then
      call shellmitc_basisfromcovariant(reference_basis, material_local_basis, material_reciprocal_basis, material_jacobian)
    endif
 
    dstrain = (/ strain(1), strain(2), 0.0d0, strain(3), strain(4), strain(5) /)
  end subroutine shellmitc_evaluateassumedstrain
 
  !--------------------------------------------------------------------
  subroutine shellmitc_materialtransform(local_basis, reciprocal_basis, transform, transform_inverse)
    real(kind=kreal), intent(in) :: local_basis(3,3), reciprocal_basis(3,3)
    real(kind=kreal), intent(out) :: transform(3,3), transform_inverse(3,3)
 
    ! This is e_hat_i . g^a, consistent with the transform in LinearElastic_Shell.
    transform = matmul(transpose(local_basis), reciprocal_basis)
    ! The first two reciprocal directions are tangent to the shell surface.
    transform(3,1:2) = 0.0d0
    transform_inverse = transform
    call calinverse(3, transform_inverse)
  end subroutine shellmitc_materialtransform
 
  pure subroutine shellmitc_strainvectortotensor(strain, tensor)
    real(kind=kreal), intent(in) :: strain(6)
    real(kind=kreal), intent(out) :: tensor(3,3)
 
    tensor = 0.0d0
    tensor(1,1) = strain(1)
    tensor(2,2) = strain(2)
    tensor(3,3) = strain(3)
    tensor(1,2) = 0.5d0*strain(4)
    tensor(2,1) = tensor(1,2)
    tensor(2,3) = 0.5d0*strain(5)
    tensor(3,2) = tensor(2,3)
    tensor(3,1) = 0.5d0*strain(6)
    tensor(1,3) = tensor(3,1)
  end subroutine shellmitc_strainvectortotensor
 
  pure subroutine shellmitc_tensortostrainvector(tensor, strain)
    real(kind=kreal), intent(in) :: tensor(3,3)
    real(kind=kreal), intent(out) :: strain(6)
 
    strain(1:3) = (/ tensor(1,1), tensor(2,2), tensor(3,3) /)
    strain(4:6) = 2.0d0*(/ tensor(1,2), tensor(2,3), tensor(3,1) /)
  end subroutine shellmitc_tensortostrainvector
 
  subroutine shellmitc_shellstraintolocal(strain, transform, local_strain)
    real(kind=kreal), intent(in) :: strain(6), transform(3,3)
    real(kind=kreal), intent(out) :: local_strain(6)
    real(kind=kreal) :: tensor(3,3), local_tensor(3,3)
 
    call shellmitc_strainvectortotensor(strain, tensor)
    local_tensor = matmul(transform, matmul(tensor, transpose(transform)))
    call shellmitc_tensortostrainvector(local_tensor, local_strain)
  end subroutine shellmitc_shellstraintolocal
 
  subroutine shellmitc_localstraintoshell(local_strain, transform_inverse, strain)
    real(kind=kreal), intent(in) :: local_strain(6), transform_inverse(3,3)
    real(kind=kreal), intent(out) :: strain(6)
    real(kind=kreal) :: tensor(3,3), local_tensor(3,3)
 
    call shellmitc_strainvectortotensor(local_strain, local_tensor)
    tensor = matmul(transform_inverse, matmul(local_tensor, transpose(transform_inverse)))
    call shellmitc_tensortostrainvector(tensor, strain)
  end subroutine shellmitc_localstraintoshell
 
  subroutine shellmitc_shellstresstolocal(stress, transform_inverse, local_stress)
    real(kind=kreal), intent(in) :: stress(6), transform_inverse(3,3)
    real(kind=kreal), intent(out) :: local_stress(6)
    real(kind=kreal) :: tensor(3,3), local_tensor(3,3)
 
    call shellstressvectortotensor(stress, tensor)
    local_tensor = matmul(transpose(transform_inverse), matmul(tensor, transform_inverse))
    call shelltensortostressvector(local_tensor, local_stress)
  end subroutine shellmitc_shellstresstolocal
 
  subroutine shellmitc_localstresstoshell(local_stress, transform, stress)
    real(kind=kreal), intent(in) :: local_stress(6), transform(3,3)
    real(kind=kreal), intent(out) :: stress(6)
    real(kind=kreal) :: tensor(3,3), local_tensor(3,3)
 
    call shellstressvectortotensor(local_stress, local_tensor)
    tensor = matmul(transpose(transform), matmul(local_tensor, transform))
    call shelltensortostressvector(tensor, stress)
  end subroutine shellmitc_localstresstoshell
 
  subroutine shellmitc_localtangenttoshell(local_tangent, transform, tangent)
    real(kind=kreal), intent(in) :: local_tangent(5,5), transform(3,3)
    real(kind=kreal), intent(out) :: tangent(5,5)
    integer :: j
    integer, parameter :: shell_component(5) = (/ 1, 2, 4, 5, 6 /)
    real(kind=kreal) :: shell_strain(6), local_strain(6)
    real(kind=kreal) :: local_stress(6), shell_stress(6)
 
    tangent = 0.0d0
    do j = 1, 5
      shell_strain = 0.0d0
      shell_strain(shell_component(j)) = 1.0d0
      call shellmitc_shellstraintolocal(shell_strain, transform, local_strain)
      local_stress = 0.0d0
      local_stress(shell_component) = matmul(local_tangent, local_strain(shell_component))
      call shellmitc_localstresstoshell(local_stress, transform, shell_stress)
      tangent(:,j) = shell_stress(shell_component)
    enddo
  end subroutine shellmitc_localtangenttoshell
 
  subroutine shellmitc_getdrillingalpha(gauss, n_layer, material_local_basis, &
      material_reciprocal_basis, alpha)
    use mmechgauss, only: tgaussstatus
    use m_elasticlinear, only: linearelastic_shell
    use mmaterial, only: shell
    implicit none
 
    type(tgaussstatus), intent(in) :: gauss
    integer, intent(in) :: n_layer
    real(kind=kreal), intent(in) :: material_local_basis(3,3), material_reciprocal_basis(3,3)
    real(kind=kreal), intent(out) :: alpha
 
    real(kind=kreal) :: elastic_tensor(3,3,3,3)
 
    call linearelastic_shell(gauss%pMaterial, shell, elastic_tensor, material_local_basis(:, shell_xi), &
      material_local_basis(:, shell_eta), material_local_basis(:, shell_zeta), &
      material_reciprocal_basis(:, shell_xi), material_reciprocal_basis(:, shell_eta), &
      material_reciprocal_basis(:, shell_zeta), alpha, n_layer)
  end subroutine shellmitc_getdrillingalpha
 
  subroutine shellmitc_getmaterialtangent(etype, nn, gauss, n_layer, material_local_basis, &
      material_reciprocal_basis, tangent, alpha, cdsys)
    use mmechgauss, only: tgaussstatus
    use m_matmatrix, only: matlmatrix_shell
    use m_elastoplastic, only: planestressj2tangent, j2_bad_state
    use mmaterial, only: shell, iselastic
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn
    type(tgaussstatus), intent(in) :: gauss
    integer, intent(in) :: n_layer
    real(kind=kreal), intent(in) :: material_local_basis(3,3), material_reciprocal_basis(3,3)
    real(kind=kreal), intent(out) :: tangent(5,5), alpha
    real(kind=kreal), intent(in), optional :: cdsys(3, 3)
 
    integer(kind=kint) :: ierr
    real(kind=kreal) :: local_stress(6), local_tangent(5,5)
    real(kind=kreal) :: transform(3,3), transform_inverse(3,3)
 
    if( iselastic(gauss%pMaterial%mtype) ) then
      call matlmatrix_shell(gauss, shell, tangent, material_local_basis(:, shell_xi), &
        material_local_basis(:, shell_eta), material_local_basis(:, shell_zeta), &
        material_reciprocal_basis(:, shell_xi), material_reciprocal_basis(:, shell_eta), &
        material_reciprocal_basis(:, shell_zeta), alpha, n_layer, cdsys)
      return
    endif
 
    if( .not.fstr_uses_plane_stress_j2_shell(etype, nn, gauss%pMaterial) ) then
      call shellmitc_abortelastoplasticunsupported(etype)
    endif
    if( .not.associated(gauss%fstatus) .or. .not.associated(gauss%istatus) ) &
      call shellmitc_abortmaterialerror(etype, j2_bad_state)
    call shellmitc_getdrillingalpha(gauss, n_layer, material_local_basis, material_reciprocal_basis, alpha)
    call shellmitc_materialtransform(material_local_basis, material_reciprocal_basis, &
      transform, transform_inverse)
    call shellmitc_shellstresstolocal(gauss%stress, transform_inverse, local_stress)
    call planestressj2tangent(gauss%pMaterial, local_stress, gauss%istatus(1), gauss%fstatus, &
      gauss%plstrain, local_tangent, ierr)
    if( ierr /= 0 ) call shellmitc_abortmaterialerror(etype, ierr)
    call shellmitc_localtangenttoshell(local_tangent, transform, tangent)
  end subroutine shellmitc_getmaterialtangent
 
  !--------------------------------------------------------------------
  subroutine shellmitc_updatestress(etype, nn, flag, update_state, gauss, n_layer, dstrain, &
      material_local_basis, material_reciprocal_basis, stress, alpha, cdsys)
    use mmechgauss
    use m_elastoplastic, only: update_planestressj2, j2_bad_state
    use mmaterial, only: updatelag, iselastoplastic
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn, flag
    integer, intent(in) :: n_layer
    logical, intent(in) :: update_state
    type(tgaussstatus), intent(inout) :: gauss
    real(kind=kreal), intent(in) :: dstrain(6)
    real(kind=kreal), intent(in) :: material_local_basis(3, 3)
    real(kind=kreal), intent(in) :: material_reciprocal_basis(3, 3)
    real(kind=kreal), intent(out) :: stress(6), alpha
    real(kind=kreal), intent(in), optional :: cdsys(3, 3)
 
    real(kind=kreal) :: d(5, 5), strain(5), stress_shell(5)
    real(kind=kreal) :: dstress(6), dstress_trace(6), trace_coeff
    if( iselastoplastic(gauss%pMaterial%mtype) ) then
      if( .not.associated(gauss%fstatus) .or. .not.associated(gauss%istatus) ) then
        call shellmitc_abortmaterialerror(etype, j2_bad_state)
      endif
      block
        ! Perfect plasticity uses the two-entry state allocated by fstr_init_gauss.
        real(kind=kreal) :: fstatus_out(2)
        real(kind=kreal) :: shell_strain_increment(6), local_strain_increment(6), local_stress_bak(6)
        real(kind=kreal) :: local_stress(6), local_strain_increment_out(6), local_tangent(5,5)
        real(kind=kreal) :: shell_strain_increment_out(6), plpotential, energy_increment
        real(kind=kreal) :: transform(3,3), transform_inverse(3,3)
        integer(kind=kint) :: ierr
        integer :: istat
 
        shell_strain_increment = dstrain-gauss%strain_bak(1:6)
        call shellmitc_materialtransform(material_local_basis, material_reciprocal_basis, &
          transform, transform_inverse)
        call shellmitc_shellstraintolocal(shell_strain_increment, transform, local_strain_increment)
        call shellmitc_shellstresstolocal(gauss%stress_bak, transform_inverse, local_stress_bak)
        call update_planestressj2(gauss%pMaterial, local_strain_increment, local_stress_bak, &
          gauss%plstrain, gauss%fstatus, local_stress, local_strain_increment_out, local_tangent, istat, &
          fstatus_out, plpotential, energy_increment, ierr)
        if( ierr /= 0 ) call shellmitc_abortmaterialerror(etype, ierr)
        call shellmitc_localstresstoshell(local_stress, transform, dstress)
        call shellmitc_getdrillingalpha(gauss, n_layer, material_local_basis, &
          material_reciprocal_basis, alpha)
        stress = dstress
        if( .not. update_state ) return
 
        call shellmitc_localstraintoshell(local_strain_increment_out, transform_inverse, &
          shell_strain_increment_out)
        gauss%strain(1:6) = gauss%strain_bak(1:6)+shell_strain_increment_out
        gauss%stress(1:6) = dstress
        gauss%istatus(1) = istat
        gauss%fstatus = fstatus_out
        gauss%plpotential = plpotential
        gauss%strain_energy = gauss%strain_energy_bak+energy_increment
        stress = gauss%stress(1:6)
      end block
      return
    endif
 
    strain = (/ dstrain(1), dstrain(2), dstrain(4), dstrain(5), dstrain(6) /)
    call shellmitc_getmaterialtangent(etype, nn, gauss, n_layer, material_local_basis, &
      material_reciprocal_basis, d, alpha, cdsys)
    stress_shell = matmul(d, strain)
    dstress = (/ stress_shell(1), stress_shell(2), 0.0d0, stress_shell(3), stress_shell(4), stress_shell(5) /)
    stress = dstress
    if( .not. update_state ) return
 
    if( flag == updatelag ) then
      trace_coeff = shellplanestresstracecoeff(gauss, n_layer)
      call shellobjectivetracestressincrement(gauss%stress_bak(1:6), dstrain, dstress_trace, trace_coeff)
      gauss%strain(1:6) = gauss%strain_bak(1:6)+dstrain
      gauss%stress(1:6) = gauss%stress_bak(1:6)+dstress_trace+dstress
      gauss%strain_energy = gauss%strain_energy_bak+dot_product(gauss%stress(1:6), dstrain)
      gauss%strain_energy = gauss%strain_energy-0.5d0*dot_product(dstress, dstrain)
    else
      gauss%strain(1:6) = dstrain
      gauss%stress(1:6) = dstress
      gauss%strain_energy = 0.5d0*dot_product(gauss%stress(1:6), gauss%strain(1:6))
    endif
    stress = gauss%stress(1:6)
  end subroutine shellmitc_updatestress
 
  !--------------------------------------------------------------------
  subroutine shellmitc_transformoutput(use_gl_strain, S, E, covariant_basis, &
      reciprocal_basis, reference_basis, current_basis, det_ref, det_cur, strain_out, stress_out)
    use m_fstr, only: opsstype, kopss_solution
    use m_utilities, only: get_principal
    implicit none
 
    logical, intent(in) :: use_gl_strain
    real(kind=kreal), intent(in) :: s(3, 3), e(3, 3)
    real(kind=kreal), intent(in) :: covariant_basis(3, 3), reciprocal_basis(3, 3)
    real(kind=kreal), intent(in) :: reference_basis(3, 3), current_basis(3, 3)
    real(kind=kreal), intent(in) :: det_ref, det_cur
    real(kind=kreal), intent(out) :: strain_out(6), stress_out(6)
 
    integer :: i
    real(kind=kreal) :: output_basis(3, 3), output_reciprocal_basis(3, 3)
    real(kind=kreal) :: reference_reciprocal_basis(3, 3), reference_local_basis(3, 3)
    real(kind=kreal) :: stress_tensor(3, 3), strain_tensor(3, 3)
    real(kind=kreal) :: stretch_b(3, 3), tensor(6), eigval(3), princ(3, 3)
    real(kind=kreal) :: logstrain(3, 3), cg_metric(3, 3)
    real(kind=kreal) :: det, jac, eig_norm
 
    output_basis = covariant_basis
    output_reciprocal_basis = reciprocal_basis
 
    if( use_gl_strain ) then
      output_basis = reference_basis
      call shellmitc_basisfromcovariant(reference_basis, reference_local_basis, reference_reciprocal_basis, det)
      output_reciprocal_basis = reference_reciprocal_basis
    endif
 
    stress_tensor = matmul(output_basis, matmul(s, transpose(output_basis)))
    strain_tensor = matmul(output_reciprocal_basis, matmul(e, transpose(output_reciprocal_basis)))
    call shelltensortostressvector(stress_tensor, stress_out)
    strain_out(1:3) = (/ strain_tensor(1, 1), strain_tensor(2, 2), strain_tensor(3, 3) /)
    ! The finite-rotation Green-Lagrange output reports the tensor shear (no
    ! engineering factor 2), matching the original ElementStress_Shell_MITC;
    ! only the small-strain/UL output uses the engineering shear.
    if( use_gl_strain ) then
      strain_out(4:6) = (/ strain_tensor(1, 2), strain_tensor(2, 3), strain_tensor(3, 1) /)
    else
      strain_out(4:6) = 2.0d0*(/ strain_tensor(1, 2), strain_tensor(2, 3), strain_tensor(3, 1) /)
    endif
 
    if( .not. (use_gl_strain .and. opsstype == kopss_solution) ) return
 
    jac = det_cur/det_ref
    if( abs(jac) <= tiny(1.0d0) ) stop "Fail to convert shell stress: detF=0"
 
    stress_tensor = matmul(current_basis, matmul(s, transpose(current_basis)))/jac
    cg_metric = matmul(transpose(reference_reciprocal_basis), reference_reciprocal_basis)
    stretch_b = matmul(current_basis, matmul(cg_metric, transpose(current_basis)))
 
    call shelltensortostressvector(stretch_b, tensor)
    call get_principal(tensor, eigval, princ)
    do i = 1, 3
      if( eigval(i) <= 0.0d0 ) stop "Fail to calc shell log strain: stretch<0"
      eigval(i) = 0.5d0*dlog(eigval(i))
      eig_norm = dsqrt(dot_product(princ(:, i), princ(:, i)))
      if( eig_norm <= 0.0d0 ) stop "Fail to calc shell log strain: direction vector=0"
      princ(:, i) = princ(:, i)/eig_norm
    end do
 
    logstrain = 0.0d0
    do i = 1, 3
      logstrain = logstrain+eigval(i)*outer_product3(princ(:, i), princ(:, i))
    end do
 
    call shelltensortostressvector(stress_tensor, stress_out)
    strain_out(1:3) = (/ logstrain(1, 1), logstrain(2, 2), logstrain(3, 3) /)
    strain_out(4:6) = 2.0d0*(/ logstrain(1, 2), logstrain(2, 3), logstrain(3, 1) /)
 
  end subroutine shellmitc_transformoutput
 
  pure subroutine shellstressvectortotensor(stress, tensor)
    real(kind=kreal), intent(in) :: stress(6)
    real(kind=kreal), intent(out) :: tensor(3, 3)
 
    tensor = 0.0d0
    tensor(1, 1) = stress(1)
    tensor(2, 2) = stress(2)
    tensor(3, 3) = stress(3)
    tensor(1, 2) = stress(4)
    tensor(2, 1) = stress(4)
    tensor(2, 3) = stress(5)
    tensor(3, 2) = stress(5)
    tensor(3, 1) = stress(6)
    tensor(1, 3) = stress(6)
  end subroutine shellstressvectortotensor
 
  pure subroutine shelltensortostressvector(tensor, stress)
    real(kind=kreal), intent(in) :: tensor(3, 3)
    real(kind=kreal), intent(out) :: stress(6)
 
    stress(1) = tensor(1, 1)
    stress(2) = tensor(2, 2)
    stress(3) = tensor(3, 3)
    stress(4) = tensor(1, 2)
    stress(5) = tensor(2, 3)
    stress(6) = tensor(3, 1)
  end subroutine shelltensortostressvector
 
  pure subroutine shellobjectivetracestressincrement(stress_old, dstrain, dstress_trace, trace_coeff)
    real(kind=kreal), intent(in) :: stress_old(6), dstrain(6)
    real(kind=kreal), intent(out) :: dstress_trace(6)
    real(kind=kreal), intent(in), optional :: trace_coeff
 
    real(kind=kreal) :: stress_tensor(3, 3), dstress_tensor(3, 3)
    real(kind=kreal) :: trd, coeff
 
    call shellstressvectortotensor(stress_old, stress_tensor)
    coeff = 1.0d0
    if (present(trace_coeff)) coeff = trace_coeff
    trd = coeff*(dstrain(1)+dstrain(2))+dstrain(3)
    dstress_tensor = -stress_tensor*trd
    call shelltensortostressvector(dstress_tensor, dstress_trace)
  end subroutine shellobjectivetracestressincrement
 
  pure subroutine shelladdulobjectivetracetangent(stress_old, ncol, B, DB, trace_coeff)
    integer(kind=kint), intent(in) :: ncol
    real(kind=kreal), intent(in) :: stress_old(6), b(5, ncol)
    real(kind=kreal), intent(inout) :: db(5, ncol)
    real(kind=kreal), intent(in), optional :: trace_coeff
 
    integer(kind=kint) :: j
    real(kind=kreal) :: dstrain_col(6), dstress_trace(6)
 
    do j = 1, ncol
      dstrain_col = 0.0d0
      dstrain_col(1:2) = b(1:2, j)
      call shellobjectivetracestressincrement(stress_old, dstrain_col, dstress_trace, trace_coeff)
      db(1, j) = db(1, j)+dstress_trace(1)
      db(2, j) = db(2, j)+dstress_trace(2)
      db(3, j) = db(3, j)+dstress_trace(4)
      db(4, j) = db(4, j)+dstress_trace(5)
      db(5, j) = db(5, j)+dstress_trace(6)
    end do
  end subroutine shelladdulobjectivetracetangent
 
  pure subroutine shelldirectorincrementalsecondderiv(director_current, director_second)
    real(kind=kreal), intent(in) :: director_current(3)
    real(kind=kreal), intent(out) :: director_second(3, 3, 3)
 
    integer :: m, n
    real(kind=kreal) :: basis_m(3), basis_n(3)
    real(kind=kreal) :: cross_n(3), cross_mn(3), cross_m(3), cross_nm(3)
 
    do n = 1, 3
      basis_n = 0.0d0
      basis_n(n) = 1.0d0
      cross_n = matmul(shellskewmatrix(basis_n), director_current)
      do m = 1, 3
        basis_m = 0.0d0
        basis_m(m) = 1.0d0
        cross_m = matmul(shellskewmatrix(basis_m), director_current)
        cross_mn = matmul(shellskewmatrix(basis_m), cross_n)
        cross_nm = matmul(shellskewmatrix(basis_n), cross_m)
        director_second(:, m, n) = 0.5d0*(cross_mn+cross_nm)
      end do
    end do
  end subroutine shelldirectorincrementalsecondderiv
 
 
  subroutine shellmitc_abortnonlinearunsupported(etype)
    integer(kind=kint), intent(in) :: etype
 
    !$omp critical
    write(*,*) '###ERROR### : Element type not supported for nonlinear static analysis'
    write(*,*) ' ic_type = ', etype
    !$omp end critical
    call hecmw_abort(hecmw_comm_get_comm())
  end subroutine shellmitc_abortnonlinearunsupported
 
  subroutine shellmitc_abortelastoplasticunsupported(etype)
    integer(kind=kint), intent(in) :: etype
 
    !$omp critical
    write(*,*) '###ERROR### : Elastoplastic shell support is limited to MITC4 J2 small deformation'
    write(*,*) '                 with perfect plasticity'
    write(*,*) ' ic_type = ', etype
    !$omp end critical
    call hecmw_abort(hecmw_comm_get_comm())
  end subroutine shellmitc_abortelastoplasticunsupported
 
  subroutine shellmitc_abortmaterialerror(etype, ierr)
    use m_elastoplastic, only: planestressj2errormessage
    integer(kind=kint), intent(in) :: etype, ierr
 
    !$omp critical
    write(*,*) '###ERROR### : ', trim(planestressj2errormessage(ierr))
    write(*,*) ' ic_type = ', etype, ', material error = ', ierr
    !$omp end critical
    call hecmw_abort(hecmw_comm_get_comm())
  end subroutine shellmitc_abortmaterialerror
 
  !--------------------------------------------------------------------
  subroutine shellmitc_integratestiffnesslayer(etype, nn, ndof, ilayer, kinematics, &
      finite_rotation, use_director_tangent, use_green_lagrange, &
      add_geometric_stiffness, ecoord, elem, shell_disp, gausses, element, &
      v1, v2, v3, director, reference_director, director_tangent, &
      director_second_tangent, nodal_kinematic_dofs, nddrill, stiff_work, cdsys)
    use mmechgauss
    use m_matmatrix
    use mmaterial, only: updatelag
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn, ndof, ilayer, kinematics
    logical, intent(in) :: finite_rotation, use_director_tangent
    logical, intent(in) :: use_green_lagrange, add_geometric_stiffness
    real(kind=kreal), intent(in) :: ecoord(3, nn), elem(3, nn), shell_disp(6, nn)
    type(tgaussstatus), intent(in) :: gausses(:)
    type(telement), intent(in), optional :: element
    real(kind=kreal), intent(in) :: v1(3, nn), v2(3, nn), v3(3, nn)
    real(kind=kreal), intent(in) :: director(3, nn), reference_director(3, nn)
    real(kind=kreal), intent(in) :: director_tangent(3, 3, nn)
    real(kind=kreal), intent(in) :: director_second_tangent(3, 3, 3, nn)
    real(kind=kreal), intent(in) :: nodal_kinematic_dofs(ndof*nn)
    real(kind=kreal), intent(in), optional :: nddrill(nn)
    real(kind=kreal), intent(inout) :: stiff_work(ndof*nn, ndof*nn)
    real(kind=kreal), intent(in), optional :: cdsys(3, 3)
 
    integer :: lx, ly, ny, isize, jsize, ishell
    integer(kind=kint) :: ierr
    real(kind=kreal) :: zeta_ly, tying_zeta, xi_lx, eta_lx
    real(kind=kreal) :: surface_weight, thickness_weight, integration_jacobian
    real(kind=kreal) :: integration_weight, layer_weight
    real(kind=kreal) :: naturalcoord(2), shapefunc(nn), shapederiv(nn, 2)
    real(kind=kreal) :: d(5, 5), b(5, ndof*nn), db(5, ndof*nn)
    real(kind=kreal) :: basis_variation(3, ndof*nn, 3)
    real(kind=kreal) :: stress_vec(5), stress_old_vec(6), alpha, trace_coeff
    real(kind=kreal) :: point_triad(3, 3), covariant_basis(3, 3)
    real(kind=kreal) :: tangent_basis(3, 2), reciprocal_basis(3, 3), local_basis(3, 3)
    real(kind=kreal) :: material_reciprocal_basis(3, 3), material_local_basis(3, 3)
    real(kind=kreal) :: director_contribution(3, nn, 3)
    real(kind=kreal) :: director_second_variation(5, 3, 3, nn)
    real(kind=kreal) :: tying_b(5, ndof*nn, 6, 3)
    real(kind=kreal) :: tying_basis_variation(3, ndof*nn, 3, 6, 3)
    real(kind=kreal) :: tying_director_second_variation(5, 3, 3, nn, 6, 3)
    real(kind=kreal) :: cv(ndof*nn), cv_disp, cv_w_second(3, 3, nn)
 
    ny = numofshellthicknessquadpoints(etype)
    layer_weight = gausses(1)%pMaterial%shell_var(ilayer)%weight
    stress_vec = 0.0d0
    director_second_variation = 0.0d0
 
    ! MITC3/MITC4 use the same midsurface tying data for every thickness point.
    ! Re-evaluating it once per layer keeps the workspace layer-local by design.
    if( etype /= fe_mitc9_shell ) then
      tying_zeta = shellmitc_tyingzeta(etype, 0.0d0)
      call shellmitc_evaluatetyingstiffnessdataatzeta(etype, nn, ndof, tying_zeta, &
        elem, shell_disp, director, director_tangent, director_second_tangent, &
        use_green_lagrange, use_director_tangent, tying_b, tying_basis_variation, tying_director_second_variation)
    endif
 
    do ly = 1, ny
      call fstr_shell_layer_quadrature_gauss(etype, gausses(1), ilayer, ly, zeta_ly, thickness_weight, ierr)
      if( ierr /= 0 ) cycle
 
      ! MITC9 tying data depend on the physical thickness coordinate.
      if( etype == fe_mitc9_shell ) then
        tying_zeta = shellmitc_tyingzeta(etype, zeta_ly)
        call shellmitc_evaluatetyingstiffnessdataatzeta(etype, nn, ndof, tying_zeta, &
          elem, shell_disp, director, director_tangent, director_second_tangent, &
          use_green_lagrange, use_director_tangent, tying_b, tying_basis_variation, tying_director_second_variation)
      endif
 
      do lx = 1, numofquadpoints(etype)
        call getquadpoint(etype, lx, naturalcoord)
        xi_lx = naturalcoord(shell_xi)
        eta_lx = naturalcoord(shell_eta)
        surface_weight = getweight(etype, lx)
        call getshapefunc(etype, naturalcoord, shapefunc)
        call getshapederiv(etype, naturalcoord, shapederiv)
 
        point_triad(:, shell_xi) = matmul(v1, shapefunc)
        point_triad(:, shell_eta) = matmul(v2, shapefunc)
        point_triad(:, shell_zeta) = matmul(v3, shapefunc)
 
        call shellmitc_preparepointkinematics(etype, nn, use_green_lagrange, ecoord, elem, &
          shell_disp(1:3, :), director, reference_director, zeta_ly, shapefunc, shapederiv, &
          covariant_basis, tangent_basis, reciprocal_basis, local_basis, &
          material_reciprocal_basis, material_local_basis, integration_jacobian, director_contribution)
 
        if( present(element) .and. (add_geometric_stiffness .or. &
            fstr_uses_plane_stress_j2_shell(etype, nn, gausses(lx)%pMaterial)) ) then
          ishell = fstr_shell_layer_gauss_index(element, lx, ilayer, ly)
        else
          ishell = 0
        endif
 
        if( ishell > 0 ) then
          call shellmitc_getmaterialtangent(etype, nn, element%shell_layer_gausses(ishell), ilayer, &
            material_local_basis, material_reciprocal_basis, d, alpha, cdsys)
        else
          call shellmitc_getmaterialtangent(etype, nn, gausses(lx), ilayer, material_local_basis, &
            material_reciprocal_basis, d, alpha, cdsys)
        endif
 
        call shellmitc_buildfirststrainvariation(nn, ndof, zeta_ly, shapefunc, shapederiv, &
          covariant_basis, tangent_basis, director_contribution, director_tangent, use_director_tangent, b, basis_variation)
        if( use_director_tangent ) then
          call shellmitc_builddirectorsecondvariation(nn, zeta_ly, shapefunc, shapederiv, &
            covariant_basis, director_second_tangent, director_second_variation)
        endif
 
        if( use_director_tangent ) then
          call shellmitc_applyassumedstrain(etype, nn, ndof, xi_lx, eta_lx, tying_b, b, &
            tying_director_second_variation, director_second_variation)
        else
          call shellmitc_applyassumedstrain(etype, nn, ndof, xi_lx, eta_lx, tying_b, b)
        endif
 
        integration_weight = surface_weight*thickness_weight*integration_jacobian
 
        if( add_geometric_stiffness ) then
          if( ishell > 0 ) then
            stress_vec = (/ element%shell_layer_gausses(ishell)%stress(1), element%shell_layer_gausses(ishell)%stress(2), &
              element%shell_layer_gausses(ishell)%stress(4), element%shell_layer_gausses(ishell)%stress(5), &
              element%shell_layer_gausses(ishell)%stress(6) /)
          else
            stress_vec = (/ gausses(lx)%stress(1), gausses(lx)%stress(2), &
              gausses(lx)%stress(4), gausses(lx)%stress(5), gausses(lx)%stress(6) /)
          endif
        endif
 
        db = matmul(d, b)
        if( kinematics == updatelag .and. fstr_is_finite_rotation_shell_element(etype, nn) ) then
          if( ishell > 0 ) then
            stress_old_vec = element%shell_layer_gausses(ishell)%stress_bak(1:6)
            trace_coeff = shellplanestresstracecoeff(element%shell_layer_gausses(ishell), ilayer)
          else
            stress_old_vec = gausses(lx)%stress_bak(1:6)
            trace_coeff = shellplanestresstracecoeff(gausses(lx), ilayer)
          endif
          call shelladdulobjectivetracetangent(stress_old_vec, ndof*nn, b, db, trace_coeff=trace_coeff)
        endif
 
        do jsize = 1, ndof*nn
          do isize = 1, ndof*nn
            stiff_work(isize, jsize) = stiff_work(isize, jsize) &
              +integration_weight*layer_weight*dot_product(b(:, isize), db(:, jsize))
          end do
        end do
 
        call shellmitc_addgeometricstiffness(etype, nn, ndof, kinematics, &
          use_director_tangent, use_green_lagrange, add_geometric_stiffness, &
          xi_lx, eta_lx, integration_weight, layer_weight, stress_vec, b, &
          basis_variation, tying_basis_variation, reciprocal_basis, material_local_basis, &
          material_reciprocal_basis, director_second_variation, director, stiff_work)
 
        call shellmitc_builddrillingsecondvariation(nn, zeta_ly, shapefunc, shapederiv, &
          director_second_tangent, material_reciprocal_basis, point_triad, use_director_tangent, cv_w_second)
        call shellmitc_builddrillingvector(nn, ndof, finite_rotation, shapefunc, &
          point_triad, material_reciprocal_basis, basis_variation, nodal_kinematic_dofs, director, nddrill, cv, cv_disp)
        call shellmitc_adddrillingstiffness(nn, ndof, finite_rotation, &
          use_director_tangent, shapefunc, cv, cv_disp, cv_w_second, nodal_kinematic_dofs, &
          director, director_tangent, alpha, integration_weight, layer_weight, nddrill, stiff_work)
      end do
    end do
  end subroutine shellmitc_integratestiffnesslayer
 
  !--------------------------------------------------------------------
  subroutine shellmitc_integrateinternalforcelayer(etype, nn, ndof, ilayer, kinematics, &
      finite_rotation, use_director_tangent, use_green_lagrange, update_state, &
      ecoord, evaluation_coords, total_nodal_state, strain_nodal_state, gausses, element, &
      v1, v2, v3, director, reference_director, director_tangent, stress_elem, &
      stress_director, stress_director_increment, nodal_kinematic_dofs, nddrill, qf_work, cdsys)
    use mmechgauss
    use mmaterial, only: updatelag
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn, ndof, ilayer, kinematics
    logical, intent(in) :: finite_rotation, use_director_tangent
    logical, intent(in) :: use_green_lagrange, update_state
    real(kind=kreal), intent(in) :: ecoord(3, nn), evaluation_coords(3, nn)
    real(kind=kreal), intent(in) :: total_nodal_state(6, nn), strain_nodal_state(6, nn)
    type(tgaussstatus), intent(in) :: gausses(:)
    type(telement), intent(inout), optional :: element
    real(kind=kreal), intent(in) :: v1(3, nn), v2(3, nn), v3(3, nn)
    real(kind=kreal), intent(in) :: director(3, nn), reference_director(3, nn)
    real(kind=kreal), intent(in) :: director_tangent(3, 3, nn)
    real(kind=kreal), intent(in) :: stress_elem(3, nn), stress_director(3, nn)
    real(kind=kreal), intent(in) :: stress_director_increment(3, nn)
    real(kind=kreal), intent(in) :: nodal_kinematic_dofs(ndof*nn)
    real(kind=kreal), intent(in), optional :: nddrill(nn)
    real(kind=kreal), intent(inout) :: qf_work(ndof*nn)
    real(kind=kreal), intent(in), optional :: cdsys(3, 3)
 
    integer :: lx, ly, ny, ishell
    integer(kind=kint) :: ierr
    logical :: store_state
    type(tgaussstatus), target :: gauss_work
    type(tgaussstatus), pointer :: gauss
    real(kind=kreal) :: zeta_ly, tying_zeta, xi, eta
    real(kind=kreal) :: surface_weight, thickness_weight, integration_jacobian
    real(kind=kreal) :: integration_weight, layer_weight
    real(kind=kreal) :: naturalcoord(2), shapefunc(nn), shapederiv(nn, 2)
    real(kind=kreal) :: b(5, ndof*nn), strain_vec(5), stress_vec(5), alpha
    real(kind=kreal) :: basis_variation(3, ndof*nn, 3), point_triad(3, 3)
    real(kind=kreal) :: covariant_basis(3, 3), tangent_basis(3, 2)
    real(kind=kreal) :: reciprocal_basis(3, 3), local_basis(3, 3)
    real(kind=kreal) :: material_reciprocal_basis(3, 3), material_local_basis(3, 3)
    real(kind=kreal) :: director_contribution(3, nn, 3)
    real(kind=kreal) :: tying_b(5, ndof*nn, 6, 3), tying_strain(5, 6, 3)
    real(kind=kreal) :: dstrain(6), stress(6), strain_out(6), stress_out(6)
    real(kind=kreal) :: strain_tensor(3, 3), stress_tensor(3, 3)
    real(kind=kreal) :: stress_covariant_basis(3, 3), stress_reciprocal_basis(3, 3)
    real(kind=kreal) :: stress_material_local_basis(3, 3)
    real(kind=kreal) :: stress_material_reciprocal_basis(3, 3)
    real(kind=kreal) :: reference_basis(3, 3), current_basis(3, 3)
    real(kind=kreal) :: reference_jacobian, current_jacobian
    real(kind=kreal) :: cv(ndof*nn), cv_disp
 
    ny = numofshellthicknessquadpoints(etype)
    layer_weight = gausses(1)%pMaterial%shell_var(ilayer)%weight
    ! MITC3/MITC4 tying is evaluated at zeta=0 once for this layer.
    if( etype /= fe_mitc9_shell ) then
      tying_zeta = shellmitc_tyingzeta(etype, 0.0d0)
      call shellmitc_evaluatetyingbatzeta(etype, nn, ndof, tying_zeta, evaluation_coords, &
        total_nodal_state, director, director_tangent, use_green_lagrange, use_director_tangent, tying_b)
    endif
 
    do ly = 1, ny
      call fstr_shell_layer_quadrature_gauss(etype, gausses(1), ilayer, ly, zeta_ly, thickness_weight, ierr)
      if( ierr /= 0 ) cycle
 
      if( etype == fe_mitc9_shell ) then
        tying_zeta = shellmitc_tyingzeta(etype, zeta_ly)
        call shellmitc_evaluatetyingbatzeta(etype, nn, ndof, tying_zeta, evaluation_coords, &
          total_nodal_state, director, director_tangent, use_green_lagrange, use_director_tangent, tying_b)
      endif
 
      if( update_state ) then
        call shellmitc_evaluatetyingstrains(etype, nn, zeta_ly, stress_elem, &
          strain_nodal_state, stress_director, stress_director_increment, use_green_lagrange, tying_strain)
      endif
 
      do lx = 1, numofquadpoints(etype)
        call getquadpoint(etype, lx, naturalcoord)
        xi = naturalcoord(shell_xi)
        eta = naturalcoord(shell_eta)
        surface_weight = getweight(etype, lx)
        call getshapefunc(etype, naturalcoord, shapefunc)
        call getshapederiv(etype, naturalcoord, shapederiv)
 
        ishell = 0
        store_state = .false.
        if( update_state .and. present(element) ) then
          ishell = fstr_shell_layer_gauss_index(element, lx, ilayer, ly)
          store_state = ishell > 0
        endif
        if( update_state .and. .not. store_state ) then
          stop "Missing shell layer Gauss state"
        endif
        if( store_state ) then
          gauss => element%shell_layer_gausses(ishell)
        else
          gauss_work = gausses(lx)
          gauss => gauss_work
        endif
 
        if( update_state ) then
          call shellmitc_evaluateassumedstrain(etype, nn, naturalcoord, zeta_ly, stress_elem, &
            strain_nodal_state, stress_director, stress_director_increment, &
            use_green_lagrange, tying_strain, dstrain, strain_tensor, stress_covariant_basis, stress_reciprocal_basis, &
            stress_material_local_basis, stress_material_reciprocal_basis, &
            reference_basis, current_basis, reference_jacobian, current_jacobian)
        endif
 
        point_triad(:, shell_xi) = matmul(v1, shapefunc)
        point_triad(:, shell_eta) = matmul(v2, shapefunc)
        point_triad(:, shell_zeta) = matmul(v3, shapefunc)
        call shellmitc_preparepointkinematics(etype, nn, use_green_lagrange, ecoord, &
          evaluation_coords, total_nodal_state(1:3, :), director, reference_director, &
          zeta_ly, shapefunc, shapederiv, covariant_basis, tangent_basis, reciprocal_basis, &
          local_basis, material_reciprocal_basis, material_local_basis, integration_jacobian, director_contribution)
 
        call shellmitc_buildfirststrainvariation(nn, ndof, zeta_ly, shapefunc, shapederiv, &
          covariant_basis, tangent_basis, director_contribution, director_tangent, use_director_tangent, b, basis_variation)
        call shellmitc_applyassumedstrain(etype, nn, ndof, xi, eta, tying_b, b)
 
        if( .not. update_state ) then
          strain_vec = matmul(b, nodal_kinematic_dofs)
          dstrain = (/ strain_vec(1), strain_vec(2), 0.0d0, strain_vec(3), strain_vec(4), strain_vec(5) /)
          stress_material_local_basis = material_local_basis
          stress_material_reciprocal_basis = material_reciprocal_basis
        endif
 
        call shellmitc_updatestress(etype, nn, kinematics, store_state, gauss, ilayer, dstrain, &
          stress_material_local_basis, stress_material_reciprocal_basis, stress, alpha, cdsys)
 
        if( store_state ) then
          if( kinematics == updatelag ) then
            gauss%strain_out(1:6) = gauss%strain(1:6)
            gauss%stress_out(1:6) = gauss%stress(1:6)
          else
            call shellstressvectortotensor(stress, stress_tensor)
            if( fstr_uses_plane_stress_j2_shell(etype, nn, gauss%pMaterial) ) then
              call shellmitc_strainvectortotensor(gauss%strain(1:6), strain_tensor)
            endif
            call shellmitc_transformoutput(use_green_lagrange, stress_tensor, strain_tensor, &
              stress_covariant_basis, stress_reciprocal_basis, reference_basis, &
              current_basis, reference_jacobian, current_jacobian, strain_out, stress_out)
            gauss%strain_out(1:6) = strain_out
            gauss%stress_out(1:6) = stress_out
          endif
        endif
 
        stress_vec = (/ stress(1), stress(2), stress(4), stress(5), stress(6) /)
        integration_weight = surface_weight*thickness_weight*integration_jacobian
        qf_work = qf_work+integration_weight*layer_weight*matmul(stress_vec, b)
 
        call shellmitc_builddrillingvector(nn, ndof, finite_rotation, shapefunc, &
          point_triad, material_reciprocal_basis, basis_variation, nodal_kinematic_dofs, director, nddrill, cv, cv_disp)
        qf_work = qf_work+integration_weight*layer_weight*alpha*cv*cv_disp
        if( store_state ) then
          if( fstr_uses_plane_stress_j2_shell(etype, nn, gauss%pMaterial) ) then
            ! Quasi-Newton sums the integrated potential at the surface Gauss points.
            if( ilayer == 1 .and. ly == 1 ) element%gausses(lx)%strain_energy = 0.0d0
            element%gausses(lx)%strain_energy = element%gausses(lx)%strain_energy &
              +integration_weight*layer_weight*(gauss%strain_energy+0.5d0*alpha*cv_disp**2)
          endif
        endif
      end do
    end do
  end subroutine shellmitc_integrateinternalforcelayer
 
  !=====================================================================
  ! Public shell element entry points
  !=====================================================================
 
  !--------------------------------------------------------------------
  subroutine shellmitc_setuporientation(cdsys_ID, coords, coordsys)
    use m_common_struct, only: set_localcoordsys, g_localcoordsys
    implicit none
 
    integer(kind=kint), intent(in), optional :: cdsys_ID
    real(kind=kreal), intent(inout), optional :: coords(3, 3)
    real(kind=kreal), allocatable, intent(out) :: coordsys(:, :)
 
    integer :: serr
 
    if( .not. present(cdsys_id) .or. .not. present(coords) ) return
    if( cdsys_id <= 0 ) return
    allocate(coordsys(3, 3))
    call set_localcoordsys(coords, g_localcoordsys(cdsys_id), coordsys, serr)
    if( serr == -1 ) stop "Fail to setup local coordinate"
  end subroutine shellmitc_setuporientation
 
  !--------------------------------------------------------------------
  subroutine stf_shell_mitc(etype, nn, ndof, ecoord, gausses, stiff, thick, mixflag, &
      nddisp, element, ndtriad, ndreftriad, ndcurtriad, nddrill, cdsys_ID, coords)
 
    use mmechgauss
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn, mixflag, ndof
    real(kind=kreal), intent(in) :: ecoord(3, nn), thick
    type(tgaussstatus), intent(in) :: gausses(:)
    real(kind=kreal), intent(out) :: stiff(:, :)
    real(kind=kreal), intent(in), optional :: nddisp(ndof, nn)
    type(telement), intent(in), optional :: element
    real(kind=kreal), intent(in), optional :: ndtriad(9, nn)
    real(kind=kreal), intent(in), optional :: ndreftriad(9, nn), ndcurtriad(9, nn)
    real(kind=kreal), intent(in), optional :: nddrill(nn)
    integer(kind=kint), intent(in), optional :: cdsys_id
    real(kind=kreal), intent(inout), optional :: coords(3, 3)    
 
    integer :: ndof_shell, nb, ilayer, nlayer
    integer(kind=kint) :: kinematics
    logical :: finite_rotation, use_director_tangent, use_green_lagrange
    logical :: add_geometric_stiffness, update_state
    real(kind=kreal), allocatable :: coordsys(:, :)
    real(kind=kreal) :: elem(3, nn), shell_disp(6, nn)
    real(kind=kreal) :: v1(3, nn), v2(3, nn), v3(3, nn)
    real(kind=kreal) :: director(3, nn), reference_director(3, nn)
    real(kind=kreal) :: director_tangent(3, 3, nn)
    real(kind=kreal) :: director_second_tangent(3, 3, 3, nn)
    real(kind=kreal) :: nodal_kinematic_dofs(ndof*nn)
    real(kind=kreal) :: stiff_work(ndof*nn, ndof*nn)
 
    shell_disp = 0.0d0
    director_second_tangent = 0.0d0
    ndof_shell = min(ndof, 6)
    if( present(nddisp) ) shell_disp(1:ndof_shell, :) = nddisp(1:ndof_shell, :)
    call shellmitc_setuporientation(cdsys_id, coords, coordsys)
 
    call shellmitc_resolveformulation(etype, nn, ndof, gausses(1), present(nddisp), &
      present(element), .true., kinematics, ndof_shell, finite_rotation, &
      use_director_tangent, use_green_lagrange, add_geometric_stiffness, update_state)
 
    call shellmitc_preparenodalkinematics(etype, nn, thick, kinematics, finite_rotation, &
      use_director_tangent, use_director_tangent, ecoord, shell_disp, ndtriad, ndreftriad, &
      ndcurtriad, elem, v1, v2, v3, director, reference_director, director_tangent, director_second_tangent)
 
    stiff = 0.0d0
    stiff_work = 0.0d0
    nodal_kinematic_dofs = 0.0d0
    do nb = 1, nn
      nodal_kinematic_dofs(ndof*(nb-1)+1:ndof*(nb-1)+ndof_shell) = shell_disp(1:ndof_shell, nb)
    end do
 
    nlayer = gausses(1)%pMaterial%totallyr
    do ilayer = 1, nlayer
      call shellmitc_integratestiffnesslayer(etype, nn, ndof, ilayer, kinematics, &
        finite_rotation, use_director_tangent, use_green_lagrange, &
        add_geometric_stiffness, ecoord, elem, shell_disp, gausses, element, &
        v1, v2, v3, director, reference_director, director_tangent, &
        director_second_tangent, nodal_kinematic_dofs, nddrill, stiff_work, coordsys)
    end do
 
    stiff(1:nn*ndof, 1:nn*ndof) = stiff_work
    call shellapplymixeddofordering(mixflag, nn, ndof, stiff=stiff)
  end subroutine stf_shell_mitc
 
  subroutine stf_shell_mitc33(etype, nn, ndof, ecoord, u, u_prev, gausses, stiff, thick, mixflag, &
      element, ndtriad, ndreftriad, ndcurtriad, nddrill, cdsys_ID, coords)
    use mmechgauss
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn, ndof, mixflag
    real(kind=kreal), intent(in) :: ecoord(3, nn), u(3, nn*2), u_prev(3, nn*2), thick
    type(tgaussstatus), intent(in) :: gausses(:)
    real(kind=kreal), intent(out) :: stiff(:, :)
    type(telement), intent(in) :: element
    real(kind=kreal), intent(in) :: ndtriad(9, nn), ndreftriad(9, nn), ndcurtriad(9, nn), nddrill(nn)
    integer(kind=kint), intent(in), optional :: cdsys_id
    real(kind=kreal), intent(inout), optional :: coords(3, 3)
 
    integer :: i, sstable(ndof*nn)
    real(kind=kreal) :: mixed_disp(ndof*nn), mixed_prev(ndof*nn)
    real(kind=kreal) :: natural_disp(ndof*nn), natural_prev(ndof*nn)
    real(kind=kreal) :: shell_disp(ndof, nn), shell_prev(ndof, nn), rotation(3)
    logical :: is_mixed
 
    do i = 1, nn*2
      mixed_disp(3*(i-1)+1:3*i) = u(1:3, i)
      mixed_prev(3*(i-1)+1:3*i) = u_prev(1:3, i)
    enddo
    call shellmixeddofmap(mixflag, nn, ndof, sstable, is_mixed)
    natural_disp = mixed_disp
    natural_prev = mixed_prev
    if( is_mixed ) then
      do i = 1, ndof*nn
        natural_disp(sstable(i)) = mixed_disp(i)
        natural_prev(sstable(i)) = mixed_prev(i)
      enddo
    endif
    do i = 1, nn
      shell_disp(:, i) = natural_disp(ndof*(i-1)+1:ndof*i)
      shell_prev(:, i) = natural_prev(ndof*(i-1)+1:ndof*i)
      call shellcomposerotationvector(shell_prev(4:6, i), &
        shell_disp(4:6, i)-shell_prev(4:6, i), rotation)
      shell_disp(4:6, i) = rotation
    enddo
 
    call stf_shell_mitc(etype, nn, ndof, ecoord, gausses, stiff, thick, mixflag, &
      nddisp=shell_disp, element=element, ndtriad=ndtriad, ndreftriad=ndreftriad, &
      ndcurtriad=ndcurtriad, nddrill=nddrill, cdsys_id=cdsys_id, coords=coords)
  end subroutine stf_shell_mitc33
 
  !--------------------------------------------------------------------
  subroutine elementstress_shell_mitc(etype, nn, ndof, ecoord, gausses, edisp, &
      strain, stress, thick, zeta, n_layer, surface_gauss_points, &
      local_strain, local_stress, local_stress_override, ndtriad, ndreftriad, ndbase_disp, surface_gauss_index, &
      cdsys_ID, coords)
    use mmechgauss
    use mmaterial, only: updatelag
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn, ndof
    integer, intent(in) :: n_layer
    real(kind=kreal), intent(in) :: ecoord(3, nn), edisp(6, nn), thick, zeta
    type(tgaussstatus), intent(in) :: gausses(:)
    real(kind=kreal), intent(out) :: strain(:, :), stress(:, :)
    logical, intent(in), optional :: surface_gauss_points
    integer(kind=kint), intent(in), optional :: surface_gauss_index
    real(kind=kreal), intent(out), optional :: local_strain(:, :), local_stress(:, :)
    real(kind=kreal), intent(in), optional :: local_stress_override(:, :)
    real(kind=kreal), intent(in), optional :: ndtriad(9, nn), ndreftriad(9, nn)
    real(kind=kreal), intent(in), optional :: ndbase_disp(6, nn)
    integer(kind=kint), intent(in), optional :: cdsys_id
    real(kind=kreal), intent(inout), optional :: coords(3, 3)    
 
    integer :: lx, ip, npoints
    integer(kind=kint) :: ierr_quad, kinematics, ndof_shell
    logical :: finite_rotation, use_director_tangent, use_green_lagrange
    logical :: add_geometric_stiffness, update_state, use_surface_gauss, use_override
    type(tgaussstatus), target :: gauss_work
    real(kind=kreal), allocatable :: coordsys(:, :)
    real(kind=kreal) :: elem(3, nn), naturalcoord(2), nncoord(nn, 2), zeta_ly
    real(kind=kreal) :: director(3, nn), director_increment(3, nn)
    real(kind=kreal) :: v1(3, nn), v2(3, nn), v3(3, nn)
    real(kind=kreal) :: a_over_2_v3(3, nn), a_over_2_v3_ref(3, nn)
    real(kind=kreal) :: a_over_2_v3_deriv(3, 3, nn), a_over_2_v3_second(3, 3, 3, nn)
    real(kind=kreal) :: tying_strain(5, 6, 3)
    real(kind=kreal) :: point_strain(6), point_stress(6), alpha
    real(kind=kreal) :: strain_tensor(3, 3), stress_tensor(3, 3)
    real(kind=kreal) :: covariant_basis(3, 3), reciprocal_basis(3, 3)
    real(kind=kreal) :: material_local_basis(3, 3), material_reciprocal_basis(3, 3)
    real(kind=kreal) :: reference_basis(3, 3), current_basis(3, 3)
    real(kind=kreal) :: reference_jacobian, current_jacobian
 
    call shellmitc_setuporientation(cdsys_id, coords, coordsys)
 
    use_surface_gauss = .false.
    if( present(surface_gauss_points) ) use_surface_gauss = surface_gauss_points
    if( present(surface_gauss_index) ) then
      if( surface_gauss_index<1 .or. surface_gauss_index>numofquadpoints(etype) ) &
        stop "Invalid shell surface Gauss point"
      use_surface_gauss = .true.
    endif
 
    ! Reproduce the same stress-evaluation configuration used by UPDATE.
    call shellmitc_resolveformulation(etype, nn, ndof, gausses(1), .true., .false., &
      .false., kinematics, ndof_shell, finite_rotation, use_director_tangent, &
      use_green_lagrange, add_geometric_stiffness, update_state)
    call getnodalnaturalcoord(etype, nncoord)
    elem = ecoord
    if( kinematics == updatelag .and. present(ndbase_disp) ) elem = elem+ndbase_disp(1:3, :)
    if( kinematics == updatelag ) elem = elem+0.5d0*edisp(1:3, :)
    call shellmitc_setupnodaldirectors(etype, nn, thick, kinematics, elem, edisp, &
      finite_rotation, use_director_tangent, .false., ndtriad, ndreftriad, &
      v1=v1, v2=v2, v3=v3, a_over_2_v3=a_over_2_v3, a_over_2_v3_ref=a_over_2_v3_ref, &
      a_over_2_v3_deriv=a_over_2_v3_deriv, a_over_2_v3_second=a_over_2_v3_second)
 
    call shellmitc_setupstresskinematics(nn, ecoord, kinematics, finite_rotation, edisp, &
      a_over_2_v3, a_over_2_v3_ref, a_over_2_v3_deriv, ndbase_disp, elem, director, director_increment)
 
    call shellmitc_evaluatetyingstrains(etype, nn, zeta, elem, edisp, director, &
      director_increment, use_green_lagrange, tying_strain)
    call fstr_shell_layer_zeta(gausses(1), n_layer, zeta, zeta_ly, ierr_quad)
    if( ierr_quad /= 0 ) stop "Invalid shell layer zeta"
 
    npoints = nn
    if( use_surface_gauss ) npoints = numofquadpoints(etype)
    if( present(surface_gauss_index) ) npoints = 1
    do lx = 1, npoints
      ip = lx
      if( present(surface_gauss_index) ) ip = surface_gauss_index
      if( use_surface_gauss ) then
        call getquadpoint(etype, ip, naturalcoord)
      else
        naturalcoord = nncoord(lx, :)
      endif
 
      call shellmitc_evaluateassumedstrain(etype, nn, naturalcoord, zeta_ly, elem, edisp, &
        director, director_increment, use_green_lagrange, tying_strain, point_strain, &
        strain_tensor, covariant_basis, reciprocal_basis, material_local_basis, &
        material_reciprocal_basis, reference_basis, current_basis, reference_jacobian, current_jacobian)
 
      use_override = .false.
      if( present(local_stress_override) ) then
        use_override = ip <= size(local_stress_override, 1) .and. size(local_stress_override, 2) >= 6
      endif
      if( use_override ) then
        point_stress = (/ local_stress_override(ip, 1), local_stress_override(ip, 2), 0.0d0, local_stress_override(ip, 4), &
          local_stress_override(ip, 5), local_stress_override(ip, 6) /)
      else
        gauss_work = gausses(ip)
        call shellmitc_updatestress(etype, nn, kinematics, update_state, gauss_work, n_layer, &
          point_strain, material_local_basis, material_reciprocal_basis, point_stress, alpha, coordsys)
      endif
 
      call shellstressvectortotensor(point_stress, stress_tensor)
      call shellmitc_transformoutput(use_green_lagrange, stress_tensor, strain_tensor, &
        covariant_basis, reciprocal_basis, reference_basis, current_basis, &
        reference_jacobian, current_jacobian, strain(lx, 1:6), stress(lx, 1:6))
      if( present(local_strain) ) local_strain(lx, 1:6) = point_strain
      if( present(local_stress) ) local_stress(lx, 1:6) = point_stress
    end do
 
  end subroutine elementstress_shell_mitc
 
  !--------------------------------------------------------------------
  subroutine dl_shell(etype, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, gausses)
    use mmechgauss
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn, ndof
    real(kind=kreal), intent(in) :: xx(*), yy(*), zz(*), rho, thick, params(*)
    integer, intent(in) :: ltype
    real(kind=kreal), intent(out) :: vect(*)
    integer(kind=kint), intent(out) :: nsize
    type(tgaussstatus), intent(in) :: gausses(:)
 
    integer :: ny, lx, ly, nb, ilayer, nlayer, jsize
    integer(kind=kint) :: ierr_quad
    real(kind=kreal) :: elem(3, nn)
    real(kind=kreal) :: v1(3, nn), v2(3, nn), v3(3, nn), director(3, nn)
    real(kind=kreal) :: naturalcoord(2), shapefunc(nn), shapederiv(nn, 2)
    real(kind=kreal) :: covariant_basis(3, 3), normal_jac(3), n(3, ndof*nn)
    real(kind=kreal) :: body_force(3), position(3), origin(3), axis(3), radial(3)
    real(kind=kreal) :: val, zeta, surface_weight, thickness_weight, integration_weight, projection
    real(kind=kreal) :: layer_weight
 
    nsize = ndof*nn
    vect(1:nsize) = 0.0d0
    elem(1, :) = xx(1:nn)
    elem(2, :) = yy(1:nn)
    elem(3, :) = zz(1:nn)
    val = params(1)
    origin = params(2:4)
    axis = params(5:7)
    ny = numofshellthicknessquadpoints(etype)
 
    if (ltype >= 10) then
      do lx = 1, numofquadpoints(etype)
        call getquadpoint(etype, lx, naturalcoord)
        surface_weight = getweight(etype, lx)
        call getshapefunc(etype, naturalcoord, shapefunc)
        call getshapederiv(etype, naturalcoord, shapederiv)
        covariant_basis(:, shell_xi:shell_eta) = matmul(elem, shapederiv)
        call cross_product(covariant_basis(:, shell_xi), covariant_basis(:, shell_eta), normal_jac)
        do nb = 1, nn
          jsize = ndof*(nb-1)
          vect(jsize+1:jsize+3) = vect(jsize+1:jsize+3) + surface_weight*shapefunc(nb)*normal_jac*val
        end do
      end do
    else
      call shellmitc_setupreferencedirectors(etype, nn, thick, elem, v1, v2, v3, director)
      nlayer = gausses(1)%pMaterial%totallyr
      do ilayer = 1, nlayer
        layer_weight = gausses(1)%pMaterial%shell_var(ilayer)%weight
        do ly = 1, ny
          call fstr_shell_layer_quadrature_gauss(etype, gausses(1), ilayer, ly, zeta, thickness_weight, ierr_quad)
          if (ierr_quad /= 0) cycle
          do lx = 1, numofquadpoints(etype)
            call getquadpoint(etype, lx, naturalcoord)
            surface_weight = getweight(etype, lx)
            call getshapefunc(etype, naturalcoord, shapefunc)
            call getshapederiv(etype, naturalcoord, shapederiv)
            call shellmitc_covariantbasis(nn, elem, director, zeta, shapefunc, shapederiv, covariant_basis)
            integration_weight = surface_weight*thickness_weight *shellmitc_covariantjacobian(covariant_basis)*layer_weight
            call shellmitc_buildinterpolationmatrix(nn, ndof, zeta, shapefunc, director, n)
 
            body_force = 0.0d0
            select case (ltype)
            case (1)
              body_force(1) = val
            case (2)
              body_force(2) = val
            case (3)
              body_force(3) = val
            case (4)
              body_force = rho*val*origin
            case (5)
              position = matmul(elem, shapefunc)
              projection = dot_product(position-origin, axis)/dot_product(axis, axis)
              radial = position-(origin+projection*axis)
              body_force = rho*val*val*radial
            end select
            vect(1:nsize) = vect(1:nsize)+integration_weight*matmul(body_force, n)
          end do
        end do
      end do
    endif
  end subroutine dl_shell
 
  !--------------------------------------------------------------------
  subroutine dl_shell_33(ic_type, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, gausses)
    use mmechgauss
    implicit none
 
    integer(kind=kint), intent(in) :: ic_type, nn, ndof
    real(kind=kreal), intent(in) :: xx(*), yy(*), zz(*), rho, thick, params(*)
    integer, intent(in) :: ltype
    real(kind=kreal), intent(out) :: vect(*)
    integer(kind=kint), intent(out) :: nsize
    type(tgaussstatus), intent(in) :: gausses(:)
 
    select case (ic_type)
    case (761)
      call dl_shell(fe_mitc3_shell, 3_kint, 6_kint, xx, yy, zz, rho, thick, ltype, params, vect, nsize, gausses)
      call shellapplymixeddofordering(2_kint, 3_kint, 6_kint, qf=vect(1:18))
    case (781)
      call dl_shell(fe_mitc4_shell, 4_kint, 6_kint, xx, yy, zz, rho, thick, ltype, params, vect, nsize, gausses)
      call shellapplymixeddofordering(1_kint, 4_kint, 6_kint, qf=vect(1:24))
    case default
      nsize = nn*ndof
      vect(1:nsize) = 0.0d0
    end select
  end subroutine dl_shell_33
 
  subroutine update_shell_mitc(etype, nn, ndof, ecoord, u, du, gausses, qf, thick, mixflag, &
      nddisp, element, ndtriad, ndreftriad, ndcurtriad, nddrill, cdsys_ID, coords)
 
    use mmechgauss
    use mmaterial, only: updatelag
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn, ndof, mixflag
    real(kind=kreal), intent(in) :: ecoord(3, nn), u(:, :), du(:, :), thick
    type(tgaussstatus), intent(in) :: gausses(:)
    real(kind=kreal), intent(out) :: qf(:)
    real(kind=kreal), intent(in), optional :: nddisp(ndof, nn)
    type(telement), intent(inout), optional :: element
    real(kind=kreal), intent(in), optional :: ndtriad(9, nn), ndreftriad(9, nn)
    real(kind=kreal), intent(in), optional :: ndcurtriad(9, nn), nddrill(nn)
    integer(kind=kint), intent(in), optional :: cdsys_id
    real(kind=kreal), intent(inout), optional :: coords(3, 3)    
 
    integer :: ndof_shell, nb, ilayer, nlayer
    integer(kind=kint) :: kinematics
    logical :: finite_rotation, use_director_tangent, use_green_lagrange
    logical :: add_geometric_stiffness, update_state
    real(kind=kreal), allocatable :: coordsys(:, :)
    real(kind=kreal) :: total_nodal_state(6, nn), strain_nodal_state(6, nn)
    real(kind=kreal) :: step_base_nodal_state(6, nn)
    real(kind=kreal) :: nodal_kinematic_dofs(ndof*nn), qf_work(ndof*nn)
    real(kind=kreal) :: evaluation_coords(3, nn), stress_elem(3, nn)
    real(kind=kreal) :: v1(3, nn), v2(3, nn), v3(3, nn)
    real(kind=kreal) :: director(3, nn), reference_director(3, nn)
    real(kind=kreal) :: director_tangent(3, 3, nn)
    real(kind=kreal) :: director_second_tangent(3, 3, 3, nn)
    real(kind=kreal) :: stress_director(3, nn), stress_director_increment(3, nn)
 
    call shellmitc_setuporientation(cdsys_id, coords, coordsys)
 
    call shellmitc_resolveformulation(etype, nn, ndof, gausses(1), .true., present(element), &
      .true., kinematics, ndof_shell, finite_rotation, use_director_tangent, &
      use_green_lagrange, add_geometric_stiffness, update_state)
 
    total_nodal_state = 0.0d0
    strain_nodal_state = 0.0d0
    step_base_nodal_state = 0.0d0
    do nb = 1, nn
      total_nodal_state(1:ndof_shell, nb) = u(1:ndof_shell, nb)+du(1:ndof_shell, nb)
      strain_nodal_state(1:ndof_shell, nb) = total_nodal_state(1:ndof_shell, nb)
      step_base_nodal_state(1:ndof_shell, nb) = u(1:ndof_shell, nb)
      if( finite_rotation ) then
        call shellcomposerotationvector(u(4:6, nb), du(4:6, nb), total_nodal_state(4:6, nb))
        strain_nodal_state(4:6, nb) = total_nodal_state(4:6, nb)
      endif
    end do
    if( present(nddisp) ) total_nodal_state(1:ndof_shell, :) = nddisp(1:ndof_shell, :)
    if( kinematics == updatelag .and. update_state ) then
      strain_nodal_state = du(1:6, 1:nn)
    endif
 
    nodal_kinematic_dofs = 0.0d0
    do nb = 1, nn
      nodal_kinematic_dofs(ndof*(nb-1)+1:ndof*(nb-1)+ndof_shell) = total_nodal_state(1:ndof_shell, nb)
    end do
    qf_work = 0.0d0
    director_second_tangent = 0.0d0
    stress_elem = 0.0d0
    stress_director = 0.0d0
    stress_director_increment = 0.0d0
 
    call shellmitc_preparenodalkinematics(etype, nn, thick, kinematics, finite_rotation, &
      use_director_tangent, .false., ecoord, total_nodal_state, ndtriad, ndreftriad, &
      ndcurtriad, evaluation_coords, v1, v2, v3, director, reference_director, director_tangent, director_second_tangent)
 
    if( update_state ) then
      call shellmitc_setupstresskinematics(nn, ecoord, kinematics, finite_rotation, &
        strain_nodal_state, director, reference_director, director_tangent, &
        step_base_nodal_state, stress_elem, stress_director, stress_director_increment)
    endif
 
    nlayer = gausses(1)%pMaterial%totallyr
    do ilayer = 1, nlayer
      call shellmitc_integrateinternalforcelayer(etype, nn, ndof, ilayer, kinematics, &
        finite_rotation, use_director_tangent, use_green_lagrange, update_state, &
        ecoord, evaluation_coords, total_nodal_state, strain_nodal_state, gausses, element, &
        v1, v2, v3, director, reference_director, director_tangent, stress_elem, &
        stress_director, stress_director_increment, nodal_kinematic_dofs, nddrill, qf_work, coordsys)
    end do
 
    qf(1:ndof*nn) = qf_work
    call shellapplymixeddofordering(mixflag, nn, ndof, qf=qf(1:ndof*nn))
  end subroutine update_shell_mitc
 
  subroutine update_shell_mitc33(etype, nn, ndof, ecoord, u, du, gausses, qf, thick, mixflag, &
      nddisp, element, ndtriad, ndreftriad, ndcurtriad, nddrill, cdsys_ID, coords)
 
    use mmechgauss
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn, ndof, mixflag
    real(kind=kreal), intent(in) :: ecoord(3, nn), u(3, nn*2), du(3, nn*2), thick
    type(tgaussstatus), intent(in) :: gausses(:)
    real(kind=kreal), intent(out) :: qf(:)
    real(kind=kreal), intent(in), optional :: nddisp(3, nn)
    type(telement), intent(inout), optional :: element
    real(kind=kreal), intent(in), optional :: ndtriad(9, nn), ndreftriad(9, nn), ndcurtriad(9, nn), nddrill(nn)
    integer(kind=kint), intent(in), optional :: cdsys_id
    real(kind=kreal), intent(inout), optional :: coords(3, 3)
 
    integer :: i, sstable(ndof*nn)
    real(kind=kreal) :: mixed_u(ndof*nn), mixed_du(ndof*nn)
    real(kind=kreal) :: natural_u(ndof*nn), natural_du(ndof*nn)
    real(kind=kreal) :: shell_u(6, nn), shell_du(6, nn), shell_nddisp(6, nn)
    logical :: is_mixed
 
    do i = 1, nn*2
      mixed_u(3*(i-1)+1:3*i) = u(1:3, i)
      mixed_du(3*(i-1)+1:3*i) = du(1:3, i)
    enddo
 
    call shellmixeddofmap(mixflag, nn, ndof, sstable, is_mixed)
    natural_u = mixed_u
    natural_du = mixed_du
    if( is_mixed ) then
      do i = 1, ndof*nn
        natural_u(sstable(i)) = mixed_u(i)
        natural_du(sstable(i)) = mixed_du(i)
      end do
    endif
 
    do i = 1, nn
      shell_u(:, i) = natural_u(ndof*(i-1)+1:ndof*i)
      shell_du(:, i) = natural_du(ndof*(i-1)+1:ndof*i)
    end do
    if( present(element) ) then
      if( present(nddisp) ) then
        shell_nddisp = shell_u+shell_du
        shell_nddisp(1:3, :) = nddisp
        call update_shell_mitc(etype, nn, ndof, ecoord, shell_u, shell_du, gausses, qf, thick, mixflag, &
          nddisp=shell_nddisp, element=element, ndtriad=ndtriad, ndreftriad=ndreftriad, &
          ndcurtriad=ndcurtriad, nddrill=nddrill, cdsys_id=cdsys_id, coords=coords)
      else
        call update_shell_mitc(etype, nn, ndof, ecoord, shell_u, shell_du, gausses, qf, thick, mixflag, &
          element=element, ndtriad=ndtriad, ndreftriad=ndreftriad, &
          ndcurtriad=ndcurtriad, nddrill=nddrill, cdsys_id=cdsys_id, coords=coords)
      endif
    else
      if( present(nddisp) ) then
        shell_nddisp = shell_u+shell_du
        shell_nddisp(1:3, :) = nddisp
        call update_shell_mitc(etype, nn, ndof, ecoord, shell_u, shell_du, gausses, qf, thick, mixflag, &
          nddisp=shell_nddisp, cdsys_id=cdsys_id, coords=coords)
      else
        call update_shell_mitc(etype, nn, ndof, ecoord, shell_u, shell_du, gausses, qf, thick, mixflag, &
          cdsys_id=cdsys_id, coords=coords)
      endif
    endif
  end subroutine update_shell_mitc33
  
 
  !--------------------------------------------------------------------
  subroutine mass_shell(etype, nn, elem, rho, thick, gausses, mass, lumped)
    use mmechgauss
    implicit none
 
    integer(kind=kint), intent(in) :: etype, nn
    real(kind=kreal), intent(in) :: elem(3, nn), rho, thick
    type(tgaussstatus), intent(in) :: gausses(:)
    real(kind=kreal), intent(out) :: mass(:, :), lumped(:)
 
    integer(kind=kint), parameter :: ndof = 6_kint
    integer :: ny, nsize, lx, ly, nb, i, ilayer, nlayer, idof
    integer(kind=kint) :: ierr_quad
    real(kind=kreal) :: v1(3, nn), v2(3, nn), v3(3, nn), director(3, nn)
    real(kind=kreal) :: naturalcoord(2), shapefunc(nn), shapederiv(nn, 2)
    real(kind=kreal) :: covariant_basis(3, 3), n(3, ndof*nn)
    real(kind=kreal) :: zeta, surface_weight, thickness_weight, integration_weight, layer_weight
    real(kind=kreal) :: totalmass, totdiag
 
    nsize = ndof*nn
    mass(1:nsize, 1:nsize) = 0.0d0
    lumped(1:nsize) = 0.0d0
    totalmass = 0.0d0
    ny = numofshellthicknessquadpoints(etype)
    call shellmitc_setupreferencedirectors(etype, nn, thick, elem, v1, v2, v3, director)
 
    nlayer = gausses(1)%pMaterial%totallyr
    do ilayer = 1, nlayer
      layer_weight = gausses(1)%pMaterial%shell_var(ilayer)%weight
      do ly = 1, ny
        call fstr_shell_layer_quadrature_gauss(etype, gausses(1), ilayer, ly, zeta, thickness_weight, ierr_quad)
        if (ierr_quad /= 0) cycle
        do lx = 1, numofquadpoints(etype)
          call getquadpoint(etype, lx, naturalcoord)
          surface_weight = getweight(etype, lx)
          call getshapefunc(etype, naturalcoord, shapefunc)
          call getshapederiv(etype, naturalcoord, shapederiv)
          call shellmitc_covariantbasis(nn, elem, director, zeta, shapefunc, shapederiv, covariant_basis)
          integration_weight = surface_weight*thickness_weight *shellmitc_covariantjacobian(covariant_basis)*layer_weight
          call shellmitc_buildinterpolationmatrix(nn, ndof, zeta, shapefunc, director, n)
          mass(1:nsize, 1:nsize) = mass(1:nsize, 1:nsize) + rho*integration_weight*matmul(transpose(n), n)
          totalmass = totalmass+rho*integration_weight
        end do
      end do
    end do
 
    totalmass = 3.0d0*totalmass
    totdiag = 0.0d0
    do nb = 1, nn
      do i = 1, 3
        idof = ndof*(nb-1)+i
        totdiag = totdiag+mass(idof, idof)
      end do
    end do
    do nb = 1, nn
      do i = 1, ndof
        idof = ndof*(nb-1)+i
        lumped(idof) = mass(idof, idof)*totalmass/totdiag
      end do
    end do
  end subroutine mass_shell
 
end module m_static_lib_shell