dpsi_square_second Module Function

    pure elemental module function dpsi_square_second(self, i, j, r) result(dpsi)
        implicit none
        class(type_square_second), intent(in) :: self
        integer(int32), intent(in) :: i
        integer(int32), intent(in) :: j
        type(type_dp_vector_3d), intent(in) :: r
        real(real64) :: dpsi

        select case (j)
        case (1)
            select case (i)
            case (1)
                dpsi = 0.25d0 * (1.0d0 - r%y) * (2.0d0 * r%x + r%y)
            case (2)
                dpsi = 0.25d0 * (1.0d0 - r%y) * (2.0d0 * r%x - r%y)
            case (3)
                dpsi = 0.25d0 * (1.0d0 + r%y) * (2.0d0 * r%x + r%y)
            case (4)
                dpsi = 0.25d0 * (1.0d0 + r%y) * (2.0d0 * r%x - r%y)
            case (5)
                dpsi = -r%x * (1.0d0 - r%y)
            case (6)
                dpsi = 0.5d0 * (1.0d0 + r%y) * (1.0d0 - r%y)
            case (7)
                dpsi = -r%x * (1.0d0 + r%y)
            case (8)
                dpsi = -0.5d0 * (1.0d0 + r%y) * (1.0d0 - r%y)
            case default
                dpsi = 0.0d0
            end select
        case (2)
            select case (i)
            case (1)
                dpsi = 0.25d0 * (1.0d0 - r%y) * (r%x + 2.0d0 * r%y)
            case (2)
                dpsi = 0.25d0 * (1.0d0 - r%y) * (-r%x + 2.0d0 * r%y)
            case (3)
                dpsi = 0.25d0 * (1.0d0 + r%y) * (r%x + 2.0d0 * r%y)
            case (4)
                dpsi = 0.25d0 * (1.0d0 + r%y) * (-r%x + 2.0d0 * r%y)
            case (5)
                dpsi = -0.5d0 * (1.0d0 + r%x) * (1.0d0 - r%x)
            case (6)
                dpsi = -(1.0d0 + r%x) * r%y
            case (7)
                dpsi = 0.5d0 * (1.0d0 + r%x) * (1.0d0 - r%x)
            case (8)
                dpsi = -(1.0d0 - r%x) * r%y
            case default
                dpsi = 0.0d0
            end select
        end select
    end function dpsi_square_second