! Copyright (c) 2022-2026 Jason Christopherson
! SPDX-License-Identifier: MIT
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The Software is provided "as is", without warranty of any kind, express or
! implied, including but not limited to the warranties of merchantability,
! fitness for a particular purpose and noninfringement.

! References:
! - J.-L. Batoz, K.-J. Bathe, L.-W. Ho, "A study of three-node triangular
!   plate bending elements," Int. J. Numer. Meth. Engng., 15, 1771-1812, 1980.
! - E. N. Dvorkin, K.-J. Bathe, "A continuum mechanics based four-node shell
!   element for general nonlinear analysis," Eng. Comput., 1, 77-88, 1984.
! - T. J. R. Hughes, F. Brezzi, "On drilling degrees of freedom," Comput.
!   Methods Appl. Mech. Engrg., 72, 105-121, 1989.
! - R. D. Cook, D. S. Malkus, M. E. Plesha, R. J. Witt, "Concepts and
!   Applications of Finite Element Analysis," 4th ed., Wiley, 2002.

module dynamics_shell_elements
    !! Flat, three-dimensional shell elements with six degrees of freedom per
    !! node.
    !!
    !! Each element is formulated in a local, planar coordinate system by
    !! superposing a membrane (plane-stress) element, a plate-bending element,
    !! and a drilling-rotation stabilization term.  The resulting local
    !! matrices are then transformed into the global coordinate system.
    !!
    !! The nodal degrees of freedom are ordered as
    !! $$ \boldsymbol{u}_i = \begin{bmatrix} u_i & v_i & w_i & \theta_{x,i} &
    !! \theta_{y,i} & \theta_{z,i} \end{bmatrix}^T, $$
    !! where the rotations follow the right-hand rule about each axis.  The
    !! element displacement vector stacks the nodal vectors in local node
    !! order.
    !!
    !! **Local Coordinate System**
    !!
    !! The local \(z\)-axis is the unit normal of the element plane computed
    !! from the node ordering (right-hand rule).  The local \(x\)-axis is the
    !! projection of the vector from node 1 to node 2 onto the element plane,
    !! and the local \(y\)-axis completes the right-handed system.  Nodes are
    !! assumed to lie in a common plane; warped quadrilaterals are projected
    !! onto their mean plane.
    !!
    !! **Kinematics**
    !!
    !! Within the local system, the displacement through the thickness is
    !! \(u = u_0 + z\theta_y\) and \(v = v_0 - z\theta_x\).  The generalized
    !! strain vector is
    !! $$ \boldsymbol{\varepsilon} = \begin{bmatrix} \varepsilon_{xx} &
    !! \varepsilon_{yy} & \gamma_{xy} & \kappa_{xx} & \kappa_{yy} &
    !! \kappa_{xy} & \gamma_{xz} & \gamma_{yz} \end{bmatrix}^T, $$
    !! with
    !! $$ \kappa_{xx} = \frac{\partial \theta_y}{\partial x}, \quad
    !! \kappa_{yy} = -\frac{\partial \theta_x}{\partial y}, \quad
    !! \kappa_{xy} = \frac{\partial \theta_y}{\partial y} -
    !! \frac{\partial \theta_x}{\partial x}, $$
    !! $$ \gamma_{xz} = \frac{\partial w}{\partial x} + \theta_y, \quad
    !! \gamma_{yz} = \frac{\partial w}{\partial y} - \theta_x. $$
    !!
    !! The corresponding stress-resultant vector is
    !! $$ \boldsymbol{\sigma} = \begin{bmatrix} N_{xx} & N_{yy} & N_{xy} &
    !! M_{xx} & M_{yy} & M_{xy} & Q_x & Q_y \end{bmatrix}^T = D
    !! \boldsymbol{\varepsilon}, $$
    !! where \(N\) are membrane forces per unit length, \(M\) are moments per
    !! unit length, and \(Q\) are transverse shear forces per unit length.
    !!
    !! **Drilling Rotation**
    !!
    !! The in-plane rotation \(\theta_z\) is not part of classical shell
    !! theory.  It is stabilized with the Hughes-Brezzi penalty
    !! $$ \Pi_d = \frac{1}{2}\int_A \alpha G t \left( \theta_z -
    !! \frac{1}{2}\left( \frac{\partial v}{\partial x} -
    !! \frac{\partial u}{\partial y} \right) \right)^2 dA, $$
    !! which vanishes for rigid-body motion and prevents a singular stiffness
    !! matrix for co-planar meshes.  The factor \(\alpha\) is the element's
    !! drilling_factor.
    !!
    !! **Mass Matrix**
    !!
    !! The consistent mass matrix is computed from the linear (triangle) or
    !! bilinear (quadrilateral) interpolation of all six nodal quantities with
    !! translational inertia \(\rho t\) and rotary inertia
    !! \(\rho t^3 / 12\).  The rotary inertia is also assigned to the drilling
    !! rotation to keep the mass matrix positive definite.
    use iso_fortran_env, only : int32, real64
    use dynamics_error_handling
    use dynamics_helper, only : cross_product
    use dynamics_structural, only : node, material, element
    implicit none
    private
    public :: shell_element
    public :: triangular_shell_element
    public :: rectangular_shell_element

! ******************************************************************************
! TYPES
! ------------------------------------------------------------------------------
    type, extends(element), abstract :: shell_element
        !! Defines a flat shell element with six degrees of freedom per node
        !! (three translations and three rotations).
        real(real64) :: thickness
            !! The shell thickness.
        real(real64) :: drilling_factor = 1.0d-3
            !! The nondimensional drilling-rotation penalty factor
            !! \(\alpha\).  The drilling stiffness is \(\alpha G t\), where
            !! \(G\) is the shear modulus.
        real(real64) :: shear_correction = 5.0d0 / 6.0d0
            !! The transverse shear correction factor.  This value is only
            !! used by elements that include transverse shear deformation.
    contains
        procedure(shell_natural_gradient), deferred, public, pass :: &
            shape_function_natural_gradient
        procedure(shell_integration_rule), deferred, public, pass :: &
            integration_rule
        procedure, public :: get_dimensionality => shl_dimensionality
        procedure, public :: get_dof_per_node => shl_dof_per_node
        procedure, public :: constitutive_matrix => shl_constitutive_matrix
        procedure, public :: jacobian => shl_jacobian
        procedure, public :: shape_function_matrix => shl_shape_function_matrix
        procedure, public :: shape_function_gradient => &
            shl_shape_function_gradient
        procedure, public :: local_frame => shl_local_frame
        procedure, public :: local_coordinates => shl_local_coordinates
        procedure, public :: rotation_matrix => shl_rotation_matrix
        procedure, public :: area => shl_area
        procedure, public :: stiffness_matrix => shl_stiffness_matrix
        procedure, public :: mass_matrix => shl_mass_matrix
        procedure, public :: external_force_vector => shl_ext_force_vector
        procedure, public :: strain => shl_strain
        procedure, public :: stress => shl_stress
    end type

! ------------------------------------------------------------------------------
    type, extends(shell_element) :: triangular_shell_element
        !! Defines a three-node flat shell element.
        !!
        !! The membrane behavior is modeled by the constant-strain triangle
        !! (CST) and the bending behavior by the Discrete Kirchhoff Triangle
        !! (DKT) of Batoz, Bathe, and Ho.  The DKT element enforces the
        !! Kirchhoff hypothesis at discrete points; therefore, transverse
        !! shear deformation is neglected and the transverse shear strain and
        !! force components are always zero.  This element is appropriate
        !! for thin shells.
        !!
        !! The natural coordinates \((\xi, \eta)\) are the area coordinates
        !! of nodes 2 and 3, so nodes 1, 2, and 3 are located at \((0,0)\),
        !! \((1,0)\), and \((0,1)\), respectively.  Nodes should be ordered
        !! counter-clockwise when viewed from the positive local \(z\)-axis;
        !! the local \(z\)-axis is defined by this ordering.
        type(node), dimension(3) :: nodes
            !! The element nodes.
    contains
        procedure, public :: get_node_count => tri_get_node_count
        procedure, public :: get_node => tri_get_node
        procedure, public :: get_node_natural_coordinates => &
            tri_get_node_natural_coordinates
        procedure, public :: evaluate_shape_function => tri_shape_function
        procedure, public :: shape_function_natural_gradient => &
            tri_shape_function_natural_gradient
        procedure, public :: integration_rule => tri_integration_rule
        procedure, public :: strain_displacement_matrix => &
            tri_strain_disp_matrix
    end type

    interface triangular_shell_element
        module procedure :: tri_init
    end interface

! ------------------------------------------------------------------------------
    type, extends(shell_element) :: rectangular_shell_element
        !! Defines a four-node flat shell element.
        !!
        !! The membrane behavior is modeled by the bilinear isoparametric
        !! quadrilateral and the bending behavior by the Mindlin-Reissner
        !! MITC4 formulation of Bathe and Dvorkin, which uses assumed
        !! transverse shear strains tied at the edge midpoints to eliminate
        !! shear locking.  Both thick and thin shells are supported.
        !! Although nominally rectangular, the isoparametric formulation is
        !! valid for any flat, convex quadrilateral.
        !!
        !! The natural coordinates \((r, s)\) span \([-1, 1]^2\), with nodes
        !! 1 through 4 located at \((-1,-1)\), \((1,-1)\), \((1,1)\), and
        !! \((-1,1)\), respectively.  Nodes should be ordered
        !! counter-clockwise when viewed from the positive local \(z\)-axis;
        !! the local \(z\)-axis is defined by this ordering.
        type(node), dimension(4) :: nodes
            !! The element nodes.
    contains
        procedure, public :: get_node_count => quad_get_node_count
        procedure, public :: get_node => quad_get_node
        procedure, public :: get_node_natural_coordinates => &
            quad_get_node_natural_coordinates
        procedure, public :: evaluate_shape_function => quad_shape_function
        procedure, public :: shape_function_natural_gradient => &
            quad_shape_function_natural_gradient
        procedure, public :: integration_rule => quad_integration_rule
        procedure, public :: strain_displacement_matrix => &
            quad_strain_disp_matrix
    end type

    interface rectangular_shell_element
        module procedure :: quad_init
    end interface

! ******************************************************************************
! INTERFACES
! ------------------------------------------------------------------------------
    interface
        pure function shell_natural_gradient(this, s) result(rst)
            !! Defines the signature of a routine returning the derivatives of
            !! the element shape functions with respect to the natural
            !! coordinates.
            import :: shell_element, real64
            class(shell_element), intent(in) :: this
                !! The shell_element object.
            real(real64), intent(in), dimension(:) :: s
                !! The natural coordinates at which to evaluate the
                !! derivatives.
            real(real64), allocatable, dimension(:,:) :: rst
                !! A 2-by-N matrix, where N is the number of nodes, whose
                !! first row contains the derivatives with respect to the
                !! first natural coordinate and whose second row contains the
                !! derivatives with respect to the second natural coordinate.
        end function

        pure subroutine shell_integration_rule(this, pts, wts)
            !! Defines the signature of a routine returning the numerical
            !! integration rule used by the element.  An integral over the
            !! element is evaluated as
            !! \(\int_A f\,dA \approx \sum_i w_i f(\boldsymbol{s}_i)
            !! \det J(\boldsymbol{s}_i)\).
            import :: shell_element, real64
            class(shell_element), intent(in) :: this
                !! The shell_element object.
            real(real64), allocatable, intent(out), dimension(:,:) :: pts
                !! A 2-by-M matrix containing the natural coordinates of the
                !! M integration points.
            real(real64), allocatable, intent(out), dimension(:) :: wts
                !! The M integration weights.
        end subroutine
    end interface

contains
! ******************************************************************************
! PRIVATE HELPERS
! ------------------------------------------------------------------------------
pure function det2(x) result(rst)
    ! Computes the determinant of a 2-by-2 matrix.
    real(real64), intent(in), dimension(2,2) :: x
    real(real64) :: rst
    rst = x(1,1) * x(2,2) - x(1,2) * x(2,1)
end function

! ------------------------------------------------------------------------------
pure function node_positions(elem) result(rst)
    ! Returns a 3-by-N matrix of the global nodal positions.
    class(shell_element), intent(in) :: elem
    real(real64), allocatable, dimension(:,:) :: rst

    integer(int32) :: i, n
    type(node) :: nd

    n = elem%get_node_count()
    allocate(rst(3, n))
    do i = 1, n
        nd = elem%get_node(i)
        rst(:,i) = [nd%x, nd%y, nd%z]
    end do
end function

! ------------------------------------------------------------------------------
pure function membrane_strain_rows(dndx) result(rst)
    ! Builds the 3-by-6N membrane strain-displacement rows from the 2-by-N
    ! shape function gradient in local coordinates.
    real(real64), intent(in), dimension(:,:) :: dndx
    real(real64), allocatable, dimension(:,:) :: rst

    integer(int32) :: i, c

    allocate(rst(3, 6 * size(dndx, 2)), source = 0.0d0)
    do i = 1, size(dndx, 2)
        c = 6 * (i - 1)
        rst(1,c+1) = dndx(1,i)
        rst(2,c+2) = dndx(2,i)
        rst(3,c+1) = dndx(2,i)
        rst(3,c+2) = dndx(1,i)
    end do
end function

! ------------------------------------------------------------------------------
pure function drilling_row(elem, s) result(rst)
    ! Builds the 6N-element operator giving theta_z - 0.5 * (dv/dx - du/dy)
    ! at the natural coordinate s.
    class(shell_element), intent(in) :: elem
    real(real64), intent(in), dimension(:) :: s
    real(real64), allocatable, dimension(:) :: rst

    integer(int32) :: i, c, n
    real(real64), allocatable, dimension(:,:) :: dndx

    n = elem%get_node_count()
    dndx = elem%shape_function_gradient(s)
    allocate(rst(6 * n), source = 0.0d0)
    do i = 1, n
        c = 6 * (i - 1)
        rst(c+1) = 0.5d0 * dndx(2,i)
        rst(c+2) = -0.5d0 * dndx(1,i)
        rst(c+6) = elem%evaluate_shape_function(i, s)
    end do
end function

! ******************************************************************************
! SHELL_ELEMENT MEMBERS
! ------------------------------------------------------------------------------
pure function shl_dimensionality(this) result(rst)
    !! Gets the dimensionality of the element.
    class(shell_element), intent(in) :: this
        !! The shell_element object.
    integer(int32) :: rst
        !! The dimensionality (always 3).
    rst = 3
end function

! ------------------------------------------------------------------------------
pure function shl_dof_per_node(this) result(rst)
    !! Gets the number of degrees of freedom per node.
    class(shell_element), intent(in) :: this
        !! The shell_element object.
    integer(int32) :: rst
        !! The number of degrees of freedom per node (always 6).
    rst = 6
end function

! ------------------------------------------------------------------------------
pure function shl_constitutive_matrix(this) result(rst)
    !! Computes the 8-by-8 constitutive matrix relating the generalized
    !! strains to the stress resultants.
    !! $$ D = \begin{bmatrix} t C & 0 & 0 \\ 0 & \frac{t^3}{12} C & 0 \\
    !! 0 & 0 & k G t I \end{bmatrix}, \quad
    !! C = \frac{E}{1 - \nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\
    !! 0 & 0 & \frac{1 - \nu}{2} \end{bmatrix}, $$
    !! where \(k\) is the shear correction factor and
    !! \(G = E / (2 (1 + \nu))\).
    class(shell_element), intent(in) :: this
        !! The shell_element object.
    real(real64), allocatable, dimension(:,:) :: rst
        !! The 8-by-8 constitutive matrix.

    ! Local Variables
    real(real64) :: e, nu, t, g, c(3,3)

    ! Initialization
    e = this%material%modulus
    nu = this%material%poissons_ratio
    t = this%thickness
    g = e / (2.0d0 * (1.0d0 + nu))

    ! Plane-stress material matrix
    c = reshape([1.0d0, nu, 0.0d0, nu, 1.0d0, 0.0d0, &
        0.0d0, 0.0d0, 0.5d0 * (1.0d0 - nu)], [3, 3])
    c = (e / (1.0d0 - nu**2)) * c

    ! Assemble the membrane, bending, and transverse shear blocks
    allocate(rst(8, 8), source = 0.0d0)
    rst(1:3,1:3) = t * c
    rst(4:6,4:6) = (t**3 / 12.0d0) * c
    rst(7,7) = this%shear_correction * g * t
    rst(8,8) = rst(7,7)
end function

! ------------------------------------------------------------------------------
pure function shl_local_frame(this) result(rst)
    !! Computes the direction cosine matrix of the element's local
    !! coordinate system.  The rows of the matrix are the local \(x\), \(y\),
    !! and \(z\) unit vectors expressed in global coordinates; therefore,
    !! \(\boldsymbol{x}_{local} = R \boldsymbol{x}_{global}\).
    class(shell_element), intent(in) :: this
        !! The shell_element object.
    real(real64), allocatable, dimension(:,:) :: rst
        !! The 3-by-3 direction cosine matrix.

    ! Local Variables
    integer(int32) :: i, n
    real(real64) :: nrm, scale, e1(3), e2(3), e3(3), d(3)
    real(real64), allocatable, dimension(:,:) :: p

    ! Nodal positions relative to node 1
    p = node_positions(this)
    n = size(p, 2)
    do i = n, 1, -1
        p(:,i) = p(:,i) - p(:,1)
    end do
    scale = maxval(norm2(p, dim = 1))
    if (scale <= 0.0d0) error stop DYN_INVALID_INPUT_ERROR

    ! Newell's method yields the polygon normal (exact for planar polygons
    ! and a best-fit normal for mildly warped quadrilaterals).
    e3 = 0.0d0
    do i = 2, n - 1
        e3 = e3 + cross_product(p(:,i), p(:,i+1))
    end do
    nrm = norm2(e3)
    if (nrm <= sqrt(epsilon(nrm)) * scale**2) &
        error stop DYN_INVALID_INPUT_ERROR
    e3 = e3 / nrm

    ! Local x-axis: the node 1 to node 2 vector projected onto the plane
    d = p(:,2) - dot_product(p(:,2), e3) * e3
    nrm = norm2(d)
    if (nrm <= sqrt(epsilon(nrm)) * scale) error stop DYN_INVALID_INPUT_ERROR
    e1 = d / nrm
    e2 = cross_product(e3, e1)

    allocate(rst(3, 3))
    rst(1,:) = e1
    rst(2,:) = e2
    rst(3,:) = e3
end function

! ------------------------------------------------------------------------------
pure function shl_local_coordinates(this) result(rst)
    !! Computes the in-plane nodal coordinates in the element's local
    !! coordinate system.  Node 1 is located at the local origin.
    class(shell_element), intent(in) :: this
        !! The shell_element object.
    real(real64), allocatable, dimension(:,:) :: rst
        !! An N-by-2 matrix, where N is the number of nodes, containing the
        !! local \(x\) and \(y\) coordinates of each node.

    ! Local Variables
    integer(int32) :: i, n
    real(real64) :: q(3)
    real(real64), allocatable, dimension(:,:) :: r, p

    ! Process
    r = this%local_frame()
    p = node_positions(this)
    n = size(p, 2)
    allocate(rst(n, 2))
    do i = 1, n
        q = p(:,i) - p(:,1)
        rst(i,1) = dot_product(r(1,:), q)
        rst(i,2) = dot_product(r(2,:), q)
    end do
end function

! ------------------------------------------------------------------------------
pure function shl_rotation_matrix(this) result(rst)
    !! Computes the transformation matrix \(T\) relating the global element
    !! displacement vector to the local element displacement vector such
    !! that \(\boldsymbol{u}_{local} = T \boldsymbol{u}_{global}\).  The
    !! matrix is block diagonal with the 3-by-3 direction cosine matrix
    !! repeated for the translations and rotations of each node.
    class(shell_element), intent(in) :: this
        !! The shell_element object.
    real(real64), allocatable, dimension(:,:) :: rst
        !! The 6N-by-6N transformation matrix, where N is the number of
        !! nodes.

    ! Local Variables
    integer(int32) :: i, n
    real(real64), allocatable, dimension(:,:) :: r

    ! Process
    r = this%local_frame()
    n = 2 * this%get_node_count()
    allocate(rst(3 * n, 3 * n), source = 0.0d0)
    do i = 1, n
        rst(3*i-2:3*i,3*i-2:3*i) = r
    end do
end function

! ------------------------------------------------------------------------------
pure function shl_jacobian(this, s) result(rst)
    !! Computes the 2-by-2 Jacobian matrix of the mapping from natural to
    !! local coordinates.
    !! $$ J = \begin{bmatrix} \frac{\partial x}{\partial s_1} &
    !! \frac{\partial y}{\partial s_1} \\ \frac{\partial x}{\partial s_2} &
    !! \frac{\partial y}{\partial s_2} \end{bmatrix} $$
    class(shell_element), intent(in) :: this
        !! The shell_element object.
    real(real64), intent(in), dimension(:) :: s
        !! The natural coordinates at which to evaluate the Jacobian.
    real(real64), allocatable, dimension(:,:) :: rst
        !! The 2-by-2 Jacobian matrix.
    rst = matmul(this%shape_function_natural_gradient(s), &
        this%local_coordinates())
end function

! ------------------------------------------------------------------------------
pure function shl_shape_function_gradient(this, s) result(rst)
    !! Computes the derivatives of the shape functions with respect to the
    !! local \(x\) and \(y\) coordinates,
    !! \(\nabla_{xy} N = J^{-1} \nabla_{s} N\).
    class(shell_element), intent(in) :: this
        !! The shell_element object.
    real(real64), intent(in), dimension(:) :: s
        !! The natural coordinates at which to evaluate the derivatives.
    real(real64), allocatable, dimension(:,:) :: rst
        !! A 2-by-N matrix, where N is the number of nodes, containing
        !! \(\partial N_i / \partial x\) in the first row and
        !! \(\partial N_i / \partial y\) in the second row.

    ! Local Variables
    real(real64) :: jac(2,2), jinv(2,2), dj

    ! Process
    jac = this%jacobian(s)
    dj = det2(jac)
    ! A non-positive determinant indicates an inverted or non-convex element
    if (dj <= 0.0d0) error stop DYN_INVALID_INPUT_ERROR
    jinv = reshape([jac(2,2), -jac(2,1), -jac(1,2), jac(1,1)], [2, 2]) / dj
    rst = matmul(jinv, this%shape_function_natural_gradient(s))
end function

! ------------------------------------------------------------------------------
pure function shl_shape_function_matrix(this, s) result(rst)
    !! Computes the 6-by-6N shape function matrix interpolating the local
    !! nodal quantities \([u, v, w, \theta_x, \theta_y, \theta_z]\) with the
    !! element's linear (triangle) or bilinear (quadrilateral) shape
    !! functions.  This interpolation is used to form the mass matrix and
    !! the external force vector.
    class(shell_element), intent(in) :: this
        !! The shell_element object.
    real(real64), intent(in), dimension(:) :: s
        !! The natural coordinates at which to evaluate the matrix.
    real(real64), allocatable, dimension(:,:) :: rst
        !! The 6-by-6N shape function matrix, where N is the number of nodes.

    ! Local Variables
    integer(int32) :: i, k, n
    real(real64) :: ni

    ! Process
    n = this%get_node_count()
    allocate(rst(6, 6 * n), source = 0.0d0)
    do i = 1, n
        ni = this%evaluate_shape_function(i, s)
        do k = 1, 6
            rst(k, 6 * (i - 1) + k) = ni
        end do
    end do
end function

! ------------------------------------------------------------------------------
pure function shl_area(this) result(rst)
    !! Computes the area of the element.
    class(shell_element), intent(in) :: this
        !! The shell_element object.
    real(real64) :: rst
        !! The element area.

    ! Local Variables
    integer(int32) :: i
    real(real64), allocatable, dimension(:) :: wts
    real(real64), allocatable, dimension(:,:) :: pts

    ! Process
    call this%integration_rule(pts, wts)
    rst = 0.0d0
    do i = 1, size(wts)
        rst = rst + wts(i) * det2(this%jacobian(pts(:,i)))
    end do
end function

! ------------------------------------------------------------------------------
pure function shl_stiffness_matrix(this, rule) result(rst)
    !! Computes the element stiffness matrix in the global coordinate system.
    !! The local stiffness matrix is
    !! $$ K_{local} = \int_A B^T D B \, dA + \int_A \alpha G t \,
    !! \boldsymbol{b}_d \boldsymbol{b}_d^T dA, $$
    !! where \(\boldsymbol{b}_d\) is the drilling-rotation operator.  The
    !! global matrix is \(K = T^T K_{local} T\).
    class(shell_element), intent(in) :: this
        !! The shell_element object.
    integer(int32), intent(in), optional :: rule
        !! The integration rule.  This argument is unused and is present for
        !! interface compatibility; the element uses the quadrature defined by
        !! its integration_rule routine.
    real(real64), allocatable, dimension(:,:) :: rst
        !! The 6N-by-6N stiffness matrix, where N is the number of nodes.

    ! Local Variables
    integer(int32) :: i, ndof
    real(real64) :: dj, g, kd
    real(real64), allocatable, dimension(:) :: wts, bd
    real(real64), allocatable, dimension(:,:) :: pts, b, d, t, kl

    ! Initialization
    call this%integration_rule(pts, wts)
    d = this%constitutive_matrix()
    ndof = 6 * this%get_node_count()
    g = this%material%modulus / (2.0d0 * (1.0d0 + &
        this%material%poissons_ratio))
    kd = this%drilling_factor * g * this%thickness
    allocate(kl(ndof, ndof), source = 0.0d0)

    ! Integrate the membrane, bending, shear, and drilling contributions
    do i = 1, size(wts)
        dj = wts(i) * det2(this%jacobian(pts(:,i)))
        b = this%strain_displacement_matrix(pts(:,i))
        kl = kl + dj * matmul(transpose(b), matmul(d, b))
        bd = drilling_row(this, pts(:,i))
        kl = kl + (dj * kd) * spread(bd, 2, ndof) * spread(bd, 1, ndof)
    end do

    ! Transform to global coordinates, removing round-off asymmetry
    t = this%rotation_matrix()
    rst = matmul(transpose(t), matmul(kl, t))
    rst = 0.5d0 * (rst + transpose(rst))
end function

! ------------------------------------------------------------------------------
pure function shl_mass_matrix(this, rule) result(rst)
    !! Computes the consistent element mass matrix in the global coordinate
    !! system.
    !! $$ M = T^T \left( \int_A \rho N^T \Lambda N \, dA \right) T, \quad
    !! \Lambda = \mathrm{diag}\left(t, t, t, \frac{t^3}{12},
    !! \frac{t^3}{12}, \frac{t^3}{12}\right) $$
    class(shell_element), intent(in) :: this
        !! The shell_element object.
    integer(int32), intent(in), optional :: rule
        !! The integration rule.  This argument is unused and is present for
        !! interface compatibility; the element uses the quadrature defined by
        !! its integration_rule routine.
    real(real64), allocatable, dimension(:,:) :: rst
        !! The 6N-by-6N mass matrix, where N is the number of nodes.

    ! Local Variables
    integer(int32) :: i, k, ndof
    real(real64) :: dj, lambda(6)
    real(real64), allocatable, dimension(:) :: wts
    real(real64), allocatable, dimension(:,:) :: pts, n, ln, t, ml

    ! Initialization
    call this%integration_rule(pts, wts)
    ndof = 6 * this%get_node_count()
    lambda(1:3) = this%thickness
    lambda(4:6) = this%thickness**3 / 12.0d0
    lambda = this%material%density * lambda
    allocate(ml(ndof, ndof), source = 0.0d0)

    ! Integrate
    do i = 1, size(wts)
        dj = wts(i) * det2(this%jacobian(pts(:,i)))
        n = this%shape_function_matrix(pts(:,i))
        ln = n
        do k = 1, 6
            ln(k,:) = lambda(k) * n(k,:)
        end do
        ml = ml + dj * matmul(transpose(n), ln)
    end do

    ! Transform to global coordinates, removing round-off asymmetry
    t = this%rotation_matrix()
    rst = matmul(transpose(t), matmul(ml, t))
    rst = 0.5d0 * (rst + transpose(rst))
end function

! ------------------------------------------------------------------------------
pure function shl_ext_force_vector(this, q, rule) result(rst)
    !! Computes the consistent nodal force vector, in the global coordinate
    !! system, resulting from a uniform surface traction.
    !! $$ \boldsymbol{f}_i = \left( \int_A N_i \, dA \right) \boldsymbol{q} $$
    !! Only the translational degrees of freedom receive load.
    class(shell_element), intent(in) :: this
        !! The shell_element object.
    real(real64), intent(in), dimension(:) :: q
        !! The 3-element surface traction vector (force per unit area)
        !! expressed in the global coordinate system.  For a pressure \(p\)
        !! acting along the local \(z\)-axis, supply
        !! \(\boldsymbol{q} = p \boldsymbol{e}_z\), where
        !! \(\boldsymbol{e}_z\) is the third row of the local_frame matrix.
    integer(int32), intent(in), optional :: rule
        !! The integration rule.  This argument is unused and is present for
        !! interface compatibility; the element uses the quadrature defined by
        !! its integration_rule routine.
    real(real64), allocatable, dimension(:) :: rst
        !! The 6N-element force vector, where N is the number of nodes.

    ! Local Variables
    integer(int32) :: i, j, c, n
    real(real64) :: dj
    real(real64), allocatable, dimension(:) :: wts
    real(real64), allocatable, dimension(:,:) :: pts

    ! Input Checking
    if (size(q) /= 3) error stop DYN_ARRAY_SIZE_ERROR

    ! Process
    call this%integration_rule(pts, wts)
    n = this%get_node_count()
    allocate(rst(6 * n), source = 0.0d0)
    do i = 1, size(wts)
        dj = wts(i) * det2(this%jacobian(pts(:,i)))
        do j = 1, n
            c = 6 * (j - 1)
            rst(c+1:c+3) = rst(c+1:c+3) + &
                dj * this%evaluate_shape_function(j, pts(:,i)) * q
        end do
    end do
end function

! ------------------------------------------------------------------------------
pure function shl_strain(this, displacement, s) result(rst)
    !! Computes the generalized strain vector, in the local coordinate
    !! system, at the specified natural coordinate.  The components are
    !! ordered as
    !! \([\varepsilon_{xx}, \varepsilon_{yy}, \gamma_{xy}, \kappa_{xx},
    !! \kappa_{yy}, \kappa_{xy}, \gamma_{xz}, \gamma_{yz}]\).
    class(shell_element), intent(in) :: this
        !! The shell_element object.
    real(real64), intent(in), dimension(:) :: displacement
        !! The 6N-element displacement vector in the global coordinate
        !! system, where N is the number of nodes.
    real(real64), intent(in), dimension(:) :: s
        !! The natural coordinates at which to evaluate the strain.
    real(real64), allocatable, dimension(:) :: rst
        !! The 8-element generalized strain vector.

    ! Local Variables
    real(real64), allocatable, dimension(:,:) :: t

    ! Process
    t = this%rotation_matrix()
    if (size(displacement) /= size(t, 2)) error stop DYN_ARRAY_SIZE_ERROR
    rst = matmul(this%strain_displacement_matrix(s), &
        matmul(t, displacement))
end function

! ------------------------------------------------------------------------------
pure function shl_stress(this, displacement, s) result(rst)
    !! Computes the stress-resultant vector, in the local coordinate system,
    !! at the specified natural coordinate.  The components are ordered as
    !! \([N_{xx}, N_{yy}, N_{xy}, M_{xx}, M_{yy}, M_{xy}, Q_x, Q_y]\).  The
    !! surface stresses may be recovered as
    !! \(\sigma = N / t \pm 6 M / t^2\).
    class(shell_element), intent(in) :: this
        !! The shell_element object.
    real(real64), intent(in), dimension(:) :: displacement
        !! The 6N-element displacement vector in the global coordinate
        !! system, where N is the number of nodes.
    real(real64), intent(in), dimension(:) :: s
        !! The natural coordinates at which to evaluate the stress resultants.
    real(real64), allocatable, dimension(:) :: rst
        !! The 8-element stress-resultant vector.
    rst = matmul(this%constitutive_matrix(), this%strain(displacement, s))
end function

! ******************************************************************************
! TRIANGULAR_SHELL_ELEMENT MEMBERS
! ------------------------------------------------------------------------------
pure function tri_init(mat, thickness, nd1, nd2, nd3) result(rst)
    !! Constructs a new [[triangular_shell_element]].
    class(material), intent(in) :: mat
        !! The material.
    real(real64), intent(in) :: thickness
        !! The shell thickness.  This value must be positive.
    class(node), intent(in) :: nd1
        !! The first node.
    class(node), intent(in) :: nd2
        !! The second node.
    class(node), intent(in) :: nd3
        !! The third node.
    type(triangular_shell_element) :: rst
        !! The new [[triangular_shell_element]].

    ! Local Variables
    real(real64), allocatable, dimension(:,:) :: frame

    ! Input Checking
    if (thickness <= 0.0d0) error stop DYN_INVALID_INPUT_ERROR

    ! Process
    rst%material = mat
    rst%thickness = thickness
    rst%nodes(1) = nd1
    rst%nodes(2) = nd2
    rst%nodes(3) = nd3

    ! Verify the geometry is not degenerate
    frame = rst%local_frame()
end function

! ------------------------------------------------------------------------------
pure function tri_get_node_count(this) result(rst)
    !! Gets the number of nodes in the element.
    class(triangular_shell_element), intent(in) :: this
        !! The triangular_shell_element object.
    integer(int32) :: rst
        !! The number of nodes (always 3).
    rst = 3
end function

! ------------------------------------------------------------------------------
pure function tri_get_node(this, i) result(rst)
    !! Gets the requested node from the element.
    class(triangular_shell_element), intent(in) :: this
        !! The triangular_shell_element object.
    integer(int32), intent(in) :: i
        !! The local index of the node to retrieve.
    type(node) :: rst
        !! The requested node.
    if (i < 1 .or. i > 3) error stop DYN_INDEX_OUT_OF_RANGE
    rst = this%nodes(i)
end function

! ------------------------------------------------------------------------------
pure function tri_get_node_natural_coordinates(this, i) result(rst)
    !! Returns the natural coordinates of the requested node.
    class(triangular_shell_element), intent(in) :: this
        !! The triangular_shell_element object.
    integer(int32), intent(in) :: i
        !! The local index of the node.
    real(real64), allocatable, dimension(:) :: rst
        !! The 2-element natural coordinate vector of the node.
    select case (i)
    case (1)
        rst = [0.0d0, 0.0d0]
    case (2)
        rst = [1.0d0, 0.0d0]
    case (3)
        rst = [0.0d0, 1.0d0]
    case default
        error stop DYN_INDEX_OUT_OF_RANGE
    end select
end function

! ------------------------------------------------------------------------------
pure function tri_shape_function(this, i, s) result(rst)
    !! Evaluates the i-th linear shape function,
    !! \(N_1 = 1 - \xi - \eta\), \(N_2 = \xi\), \(N_3 = \eta\).
    class(triangular_shell_element), intent(in) :: this
        !! The triangular_shell_element object.
    integer(int32), intent(in) :: i
        !! The index of the shape function to evaluate.
    real(real64), intent(in), dimension(:) :: s
        !! The natural coordinates \((\xi, \eta)\).
    real(real64) :: rst
        !! The value of the i-th shape function.
    select case (i)
    case (1)
        rst = 1.0d0 - s(1) - s(2)
    case (2)
        rst = s(1)
    case (3)
        rst = s(2)
    case default
        rst = 0.0d0
    end select
end function

! ------------------------------------------------------------------------------
pure function tri_shape_function_natural_gradient(this, s) result(rst)
    !! Computes the derivatives of the linear shape functions with respect
    !! to the natural coordinates.
    class(triangular_shell_element), intent(in) :: this
        !! The triangular_shell_element object.
    real(real64), intent(in), dimension(:) :: s
        !! The natural coordinates \((\xi, \eta)\).
    real(real64), allocatable, dimension(:,:) :: rst
        !! The 2-by-3 matrix of derivatives.
    rst = reshape([-1.0d0, -1.0d0, 1.0d0, 0.0d0, 0.0d0, 1.0d0], [2, 3])
end function

! ------------------------------------------------------------------------------
pure subroutine tri_integration_rule(this, pts, wts)
    !! Returns the three-point interior quadrature rule for triangles, which
    !! integrates quadratic polynomials exactly.
    class(triangular_shell_element), intent(in) :: this
        !! The triangular_shell_element object.
    real(real64), allocatable, intent(out), dimension(:,:) :: pts
        !! The 2-by-3 matrix of integration point natural coordinates.
    real(real64), allocatable, intent(out), dimension(:) :: wts
        !! The 3 integration weights.

    real(real64), parameter :: a = 1.0d0 / 6.0d0
    real(real64), parameter :: b = 2.0d0 / 3.0d0

    pts = reshape([a, a, b, a, a, b], [2, 3])
    wts = [a, a, a]
end subroutine

! ------------------------------------------------------------------------------
pure function dkt_curvature_matrix(xy, s) result(rst)
    ! Computes the 3-by-9 DKT curvature-displacement matrix, relating the
    ! nodal [w, theta_x, theta_y] values to [kxx, kyy, kxy], following Batoz,
    ! Bathe, and Ho (1980).  The DKT rotations are beta_x = theta_y and
    ! beta_y = -theta_x.
    real(real64), intent(in), dimension(3,2) :: xy
        ! The local nodal coordinates.
    real(real64), intent(in), dimension(:) :: s
        ! The natural coordinates (xi, eta).
    real(real64) :: rst(3,9)

    ! Edge k runs from node ii(k) to node jj(k); edges 1, 2, and 3 correspond
    ! to the midside nodes 4 (2-3), 5 (3-1), and 6 (1-2) of Batoz et al.
    integer(int32), parameter :: ii(3) = [2, 3, 1]
    integer(int32), parameter :: jj(3) = [3, 1, 2]

    ! Local Variables
    integer(int32) :: k, m
    real(real64) :: xi, eta, xij, yij, l2, x31, x12, y31, y12, twoa
    real(real64) :: a(3), b(3), c(3), d(3), e(3), dn(2,6), hx(2,9), hy(2,9)

    ! Edge geometry coefficients
    do k = 1, 3
        xij = xy(ii(k),1) - xy(jj(k),1)
        yij = xy(ii(k),2) - xy(jj(k),2)
        l2 = xij**2 + yij**2
        a(k) = -xij / l2
        b(k) = 0.75d0 * xij * yij / l2
        c(k) = (0.25d0 * xij**2 - 0.5d0 * yij**2) / l2
        d(k) = -yij / l2
        e(k) = (0.25d0 * yij**2 - 0.5d0 * xij**2) / l2
    end do

    ! Derivatives of the 6-node quadratic shape functions with respect to xi
    ! (row 1) and eta (row 2)
    xi = s(1)
    eta = s(2)
    dn(1,:) = [4.0d0 * (xi + eta) - 3.0d0, 4.0d0 * xi - 1.0d0, 0.0d0, &
        4.0d0 * eta, -4.0d0 * eta, 4.0d0 * (1.0d0 - 2.0d0 * xi - eta)]
    dn(2,:) = [4.0d0 * (xi + eta) - 3.0d0, 0.0d0, 4.0d0 * eta - 1.0d0, &
        4.0d0 * xi, 4.0d0 * (1.0d0 - xi - 2.0d0 * eta), -4.0d0 * xi]

    ! Derivatives of the beta_x (Hx) and beta_y (Hy) interpolation functions
    do m = 1, 2
        hx(m,1) = 1.5d0 * (a(3) * dn(m,6) - a(2) * dn(m,5))
        hx(m,2) = b(2) * dn(m,5) + b(3) * dn(m,6)
        hx(m,3) = dn(m,1) - c(2) * dn(m,5) - c(3) * dn(m,6)
        hx(m,4) = 1.5d0 * (a(1) * dn(m,4) - a(3) * dn(m,6))
        hx(m,5) = b(3) * dn(m,6) + b(1) * dn(m,4)
        hx(m,6) = dn(m,2) - c(3) * dn(m,6) - c(1) * dn(m,4)
        hx(m,7) = 1.5d0 * (a(2) * dn(m,5) - a(1) * dn(m,4))
        hx(m,8) = b(1) * dn(m,4) + b(2) * dn(m,5)
        hx(m,9) = dn(m,3) - c(1) * dn(m,4) - c(2) * dn(m,5)

        hy(m,1) = 1.5d0 * (d(3) * dn(m,6) - d(2) * dn(m,5))
        hy(m,2) = -dn(m,1) + e(2) * dn(m,5) + e(3) * dn(m,6)
        hy(m,3) = -hx(m,2)
        hy(m,4) = 1.5d0 * (d(1) * dn(m,4) - d(3) * dn(m,6))
        hy(m,5) = -dn(m,2) + e(3) * dn(m,6) + e(1) * dn(m,4)
        hy(m,6) = -hx(m,5)
        hy(m,7) = 1.5d0 * (d(2) * dn(m,5) - d(1) * dn(m,4))
        hy(m,8) = -dn(m,3) + e(1) * dn(m,4) + e(2) * dn(m,5)
        hy(m,9) = -hx(m,8)
    end do

    ! Map natural derivatives to local x-y derivatives
    x31 = xy(3,1) - xy(1,1)
    x12 = xy(1,1) - xy(2,1)
    y31 = xy(3,2) - xy(1,2)
    y12 = xy(1,2) - xy(2,2)
    twoa = x31 * y12 - x12 * y31
    rst(1,:) = (y31 * hx(1,:) + y12 * hx(2,:)) / twoa
    rst(2,:) = (-x31 * hy(1,:) - x12 * hy(2,:)) / twoa
    rst(3,:) = (-x31 * hx(1,:) - x12 * hx(2,:) + y31 * hy(1,:) + &
        y12 * hy(2,:)) / twoa
end function

! ------------------------------------------------------------------------------
pure function tri_strain_disp_matrix(this, s) result(rst)
    !! Computes the 8-by-18 generalized strain-displacement matrix in the
    !! local coordinate system.  Rows 1-3 contain the constant-strain
    !! membrane terms, rows 4-6 the DKT curvature terms, and rows 7-8 (the
    !! transverse shear strains) are zero as the DKT formulation neglects
    !! transverse shear deformation.
    class(triangular_shell_element), intent(in) :: this
        !! The triangular_shell_element object.
    real(real64), intent(in), dimension(:) :: s
        !! The natural coordinates \((\xi, \eta)\).
    real(real64), allocatable, dimension(:,:) :: rst
        !! The 8-by-18 strain-displacement matrix.

    ! Local Variables
    integer(int32) :: i, m
    real(real64) :: bb(3,9)
    real(real64), allocatable, dimension(:,:) :: xy

    ! Membrane terms
    allocate(rst(8, 18), source = 0.0d0)
    rst(1:3,:) = membrane_strain_rows(this%shape_function_gradient(s))

    ! Bending terms: DKT local DOF (w, theta_x, theta_y) of node i map onto
    ! element DOF 3, 4, and 5 of the node.
    xy = this%local_coordinates()
    bb = dkt_curvature_matrix(xy, s)
    do i = 1, 3
        do m = 1, 3
            rst(4:6, 6 * (i - 1) + 2 + m) = bb(:, 3 * (i - 1) + m)
        end do
    end do
end function

! ******************************************************************************
! RECTANGULAR_SHELL_ELEMENT MEMBERS
! ------------------------------------------------------------------------------
pure function quad_init(mat, thickness, nd1, nd2, nd3, nd4) result(rst)
    !! Constructs a new [[rectangular_shell_element]].
    class(material), intent(in) :: mat
        !! The material.
    real(real64), intent(in) :: thickness
        !! The shell thickness.  This value must be positive.
    class(node), intent(in) :: nd1
        !! The first node.
    class(node), intent(in) :: nd2
        !! The second node.
    class(node), intent(in) :: nd3
        !! The third node.
    class(node), intent(in) :: nd4
        !! The fourth node.
    type(rectangular_shell_element) :: rst
        !! The new [[rectangular_shell_element]].

    ! Local Variables
    real(real64), allocatable, dimension(:,:) :: frame

    ! Input Checking
    if (thickness <= 0.0d0) error stop DYN_INVALID_INPUT_ERROR

    ! Process
    rst%material = mat
    rst%thickness = thickness
    rst%nodes(1) = nd1
    rst%nodes(2) = nd2
    rst%nodes(3) = nd3
    rst%nodes(4) = nd4

    ! Verify the geometry is not degenerate
    frame = rst%local_frame()
end function

! ------------------------------------------------------------------------------
pure function quad_get_node_count(this) result(rst)
    !! Gets the number of nodes in the element.
    class(rectangular_shell_element), intent(in) :: this
        !! The rectangular_shell_element object.
    integer(int32) :: rst
        !! The number of nodes (always 4).
    rst = 4
end function

! ------------------------------------------------------------------------------
pure function quad_get_node(this, i) result(rst)
    !! Gets the requested node from the element.
    class(rectangular_shell_element), intent(in) :: this
        !! The rectangular_shell_element object.
    integer(int32), intent(in) :: i
        !! The local index of the node to retrieve.
    type(node) :: rst
        !! The requested node.
    if (i < 1 .or. i > 4) error stop DYN_INDEX_OUT_OF_RANGE
    rst = this%nodes(i)
end function

! ------------------------------------------------------------------------------
pure function quad_get_node_natural_coordinates(this, i) result(rst)
    !! Returns the natural coordinates of the requested node.
    class(rectangular_shell_element), intent(in) :: this
        !! The rectangular_shell_element object.
    integer(int32), intent(in) :: i
        !! The local index of the node.
    real(real64), allocatable, dimension(:) :: rst
        !! The 2-element natural coordinate vector of the node.

    real(real64), parameter :: rn(4) = [-1.0d0, 1.0d0, 1.0d0, -1.0d0]
    real(real64), parameter :: sn(4) = [-1.0d0, -1.0d0, 1.0d0, 1.0d0]

    if (i < 1 .or. i > 4) error stop DYN_INDEX_OUT_OF_RANGE
    rst = [rn(i), sn(i)]
end function

! ------------------------------------------------------------------------------
pure function quad_shape_function(this, i, s) result(rst)
    !! Evaluates the i-th bilinear shape function,
    !! \(N_i = \frac{1}{4}(1 + r r_i)(1 + s s_i)\).
    class(rectangular_shell_element), intent(in) :: this
        !! The rectangular_shell_element object.
    integer(int32), intent(in) :: i
        !! The index of the shape function to evaluate.
    real(real64), intent(in), dimension(:) :: s
        !! The natural coordinates \((r, s)\).
    real(real64) :: rst
        !! The value of the i-th shape function.

    real(real64), parameter :: rn(4) = [-1.0d0, 1.0d0, 1.0d0, -1.0d0]
    real(real64), parameter :: sn(4) = [-1.0d0, -1.0d0, 1.0d0, 1.0d0]

    if (i < 1 .or. i > 4) then
        rst = 0.0d0
    else
        rst = 0.25d0 * (1.0d0 + s(1) * rn(i)) * (1.0d0 + s(2) * sn(i))
    end if
end function

! ------------------------------------------------------------------------------
pure function quad_shape_function_natural_gradient(this, s) result(rst)
    !! Computes the derivatives of the bilinear shape functions with respect
    !! to the natural coordinates.
    class(rectangular_shell_element), intent(in) :: this
        !! The rectangular_shell_element object.
    real(real64), intent(in), dimension(:) :: s
        !! The natural coordinates \((r, s)\).
    real(real64), allocatable, dimension(:,:) :: rst
        !! The 2-by-4 matrix of derivatives.

    real(real64), parameter :: rn(4) = [-1.0d0, 1.0d0, 1.0d0, -1.0d0]
    real(real64), parameter :: sn(4) = [-1.0d0, -1.0d0, 1.0d0, 1.0d0]

    allocate(rst(2, 4))
    rst(1,:) = 0.25d0 * rn * (1.0d0 + s(2) * sn)
    rst(2,:) = 0.25d0 * sn * (1.0d0 + s(1) * rn)
end function

! ------------------------------------------------------------------------------
pure subroutine quad_integration_rule(this, pts, wts)
    !! Returns the 2-by-2 Gauss quadrature rule.
    class(rectangular_shell_element), intent(in) :: this
        !! The rectangular_shell_element object.
    real(real64), allocatable, intent(out), dimension(:,:) :: pts
        !! The 2-by-4 matrix of integration point natural coordinates.
    real(real64), allocatable, intent(out), dimension(:) :: wts
        !! The 4 integration weights.

    real(real64) :: g

    g = 1.0d0 / sqrt(3.0d0)
    pts = reshape([-g, -g, g, -g, g, g, -g, g], [2, 4])
    wts = [1.0d0, 1.0d0, 1.0d0, 1.0d0]
end subroutine

! ------------------------------------------------------------------------------
pure function mitc4_covariant_shear(elem, r, s, dir) result(rst)
    ! Computes the 24-element operator giving the covariant transverse shear
    ! strain at (r, s):
    !   dir = 1: gamma_r = dw/dr + theta_y * dx/dr - theta_x * dy/dr
    !   dir = 2: gamma_s = dw/ds + theta_y * dx/ds - theta_x * dy/ds
    class(rectangular_shell_element), intent(in) :: elem
    real(real64), intent(in) :: r, s
    integer(int32), intent(in) :: dir
    real(real64) :: rst(24)

    integer(int32) :: i, c
    real(real64) :: ni, jac(2,2)
    real(real64), allocatable, dimension(:,:) :: dnds

    dnds = elem%shape_function_natural_gradient([r, s])
    jac = elem%jacobian([r, s])
    rst = 0.0d0
    do i = 1, 4
        c = 6 * (i - 1)
        ni = elem%evaluate_shape_function(i, [r, s])
        rst(c+3) = dnds(dir,i)
        rst(c+4) = -ni * jac(dir,2)
        rst(c+5) = ni * jac(dir,1)
    end do
end function

! ------------------------------------------------------------------------------
pure function quad_strain_disp_matrix(this, s) result(rst)
    !! Computes the 8-by-24 generalized strain-displacement matrix in the
    !! local coordinate system.  Rows 1-3 contain the bilinear membrane
    !! terms, rows 4-6 the Mindlin-Reissner curvature terms, and rows 7-8 the
    !! MITC4 assumed transverse shear strains.
    !!
    !! The covariant transverse shear strains are sampled at the edge
    !! midpoints \(A = (0, 1)\), \(B = (-1, 0)\), \(C = (0, -1)\), and
    !! \(D = (1, 0)\) and interpolated as
    !! $$ \tilde{\gamma}_r = \frac{1}{2}(1 + s)\gamma_r^A +
    !! \frac{1}{2}(1 - s)\gamma_r^C, \quad
    !! \tilde{\gamma}_s = \frac{1}{2}(1 + r)\gamma_s^D +
    !! \frac{1}{2}(1 - r)\gamma_s^B, $$
    !! and then transformed to local Cartesian components with
    !! \([\gamma_{xz}, \gamma_{yz}]^T = J^{-1}[\tilde{\gamma}_r,
    !! \tilde{\gamma}_s]^T\).
    class(rectangular_shell_element), intent(in) :: this
        !! The rectangular_shell_element object.
    real(real64), intent(in), dimension(:) :: s
        !! The natural coordinates \((r, s)\).
    real(real64), allocatable, dimension(:,:) :: rst
        !! The 8-by-24 strain-displacement matrix.

    ! Local Variables
    integer(int32) :: i, c
    real(real64) :: jac(2,2), jinv(2,2), gr(24), gs(24)
    real(real64), allocatable, dimension(:,:) :: dndx

    ! Membrane terms
    allocate(rst(8, 24), source = 0.0d0)
    dndx = this%shape_function_gradient(s)
    rst(1:3,:) = membrane_strain_rows(dndx)

    ! Bending terms
    do i = 1, 4
        c = 6 * (i - 1)
        rst(4,c+5) = dndx(1,i)
        rst(5,c+4) = -dndx(2,i)
        rst(6,c+4) = -dndx(1,i)
        rst(6,c+5) = dndx(2,i)
    end do

    ! Assumed covariant transverse shear strains
    gr = 0.5d0 * (1.0d0 + s(2)) * &
        mitc4_covariant_shear(this, 0.0d0, 1.0d0, 1) + &
        0.5d0 * (1.0d0 - s(2)) * &
        mitc4_covariant_shear(this, 0.0d0, -1.0d0, 1)
    gs = 0.5d0 * (1.0d0 + s(1)) * &
        mitc4_covariant_shear(this, 1.0d0, 0.0d0, 2) + &
        0.5d0 * (1.0d0 - s(1)) * &
        mitc4_covariant_shear(this, -1.0d0, 0.0d0, 2)

    ! Transform to local Cartesian components
    jac = this%jacobian(s)
    jinv = reshape([jac(2,2), -jac(2,1), -jac(1,2), jac(1,1)], [2, 2]) / &
        det2(jac)
    rst(7,:) = jinv(1,1) * gr + jinv(1,2) * gs
    rst(8,:) = jinv(2,1) * gr + jinv(2,2) * gs
end function

! ------------------------------------------------------------------------------
end module