"""
    UnstructuredQuadMesh{RealT<:Real} <: AbstractMesh{2}

An unstructured (possibly curved) quadrilateral mesh.

    UnstructuredQuadMesh(filename, periodic; RealT=Float64)

All mesh information, neighbour coupling, and boundary curve information is read in
from a mesh file `filename`.

!!! warning "Experimental code"
    This mesh type is experimental and can change any time.
"""
struct UnstructuredQuadMesh{RealT<:Real, CurvedSurfaceT<:CurvedSurface{RealT}} <: AbstractMesh{2}
  filename             ::String
  n_corners            ::Int
  n_surfaces           ::Int # total number of surfaces
  n_interfaces         ::Int # number of interior surfaces
  n_boundaries         ::Int # number of surfaces on the physical boundary
  n_elements           ::Int
  polydeg              ::Int
  corners              ::Array{RealT, 2}  # [ndims, n_corners]
  neighbour_information::Array{Int, 2}  # [neighbour node/element/edge ids, n_surfaces]
  boundary_names       ::Array{String, 2} # [local sides, n_elements]
  periodicity          ::Bool
  element_node_ids     ::Array{Int, 2} # [node ids, n_elements]
  element_is_curved    ::Vector{Bool}
  surface_curves       ::Array{CurvedSurfaceT, 2} # [local sides, n_elements]
end


# constructor for an unstructured mesh read in from a file
# TODO: this mesh file parsing and construction of the mesh skeleton can likely be improved in terms
#       of performance
function UnstructuredQuadMesh(filename, periodic; RealT=Float64)

  # readin all the information from the mesh file into a string array
  file_lines = readlines(open(filename))

  # readin the number of nodes, number of interfaces, number of elements and local polynomial degree
  current_line = split(file_lines[1])
  n_corners    = parse(Int, current_line[1])
  n_surfaces   = parse(Int, current_line[2])
  n_elements   = parse(Int, current_line[3])
  mesh_polydeg = parse(Int, current_line[4])

  mesh_nnodes = mesh_polydeg + 1

  # The types of structs used in the following depend on information read from
  # the mesh file. Thus, this cannot be type stable at all. Hence, we allocate
  # the memory now and introduce a function barrier before continuing to read
  # data from the file.
  corner_nodes      = Array{RealT}(undef, (2, n_corners))
  interface_info    = Array{Int}(undef, (6, n_surfaces))
  element_node_ids  = Array{Int}(undef, (4, n_elements))
  curved_check      = Vector{Int}(undef, 4)
  cornerNodeVals    = Array{RealT}(undef, (4, 2))
  tempNodes         = Array{RealT}(undef, (4, 2))
  curve_vals        = Array{RealT}(undef, (mesh_nnodes, 2))
  element_is_curved = Array{Bool}(undef, n_elements)
  CurvedSurfaceT    = CurvedSurface{RealT}
  surface_curves    = Array{CurvedSurfaceT}(undef, (4, n_elements))
  bndy_names        = Array{String}(undef, (4, n_elements))

  # create the Chebyshev-Gauss-Lobatto nodes used to represent any curved boundaries that are
  # required to construct the sides
  cheby_nodes_, _ = chebyshev_gauss_lobatto_nodes_weights(mesh_nnodes)
  bary_weights_   = barycentric_weights(cheby_nodes_)
  cheby_nodes  = SVector{mesh_nnodes}(cheby_nodes_)
  bary_weights = SVector{mesh_nnodes}(bary_weights_)

  arrays = (; corner_nodes, interface_info, element_node_ids, curved_check,
              cornerNodeVals, tempNodes, curve_vals,
              element_is_curved, surface_curves, bndy_names)
  counters = (; n_corners, n_surfaces, n_elements)

  n_boundaries = parse_mesh_file!(arrays, RealT, CurvedSurfaceT, file_lines, counters, cheby_nodes, bary_weights)

  # get the number of internal interfaces in the mesh
  if periodic
    n_interfaces = n_surfaces
  else
    n_interfaces = n_surfaces - n_boundaries
  end

  return UnstructuredQuadMesh{RealT, CurvedSurfaceT}(
    filename, n_corners, n_surfaces, n_interfaces, n_boundaries,
    n_elements, mesh_polydeg, corner_nodes,
    interface_info, bndy_names, periodic,
    element_node_ids, element_is_curved, surface_curves)
end

function parse_mesh_file!(arrays, RealT, CurvedSurfaceT, file_lines, counters, cheby_nodes, bary_weights)
  @unpack ( corner_nodes, interface_info, element_node_ids, curved_check,
            cornerNodeVals, tempNodes, curve_vals,
            element_is_curved, surface_curves, bndy_names ) = arrays
  @unpack n_corners, n_surfaces, n_elements = counters
  mesh_nnodes = length(cheby_nodes)

  # counter to step through the mesh file line by line
  file_idx = 2

  # readin an store the nodes that dictate the corners of the elements needed to construct the
  # element geometry terms
  for j in 1:n_corners
    current_line       = split(file_lines[file_idx])
    corner_nodes[1, j] = parse(RealT, current_line[1])
    corner_nodes[2, j] = parse(RealT, current_line[2])
    file_idx += 1
  end

  # readin an store the nodes that dictate the interfaces, neighbour data, and orientations contains
  # the following:
  #    interface_info[1] = start node ID
  #    interface_info[2] = end node ID
  #    interface_info[3] = ID of the primary element
  #    interface_info[4] = ID of the secondary element (if 0 then it is a physical boundary)
  #    interface_info[5] = local side ID on the primary element
  #    interface_info[6] = local side ID on the secondary element
  # container to for the interface neighbour information and connectivity
  n_boundaries = 0
  for j in 1:n_surfaces
    current_line   = split(file_lines[file_idx])
    interface_info[1, j] = parse(Int, current_line[1])
    interface_info[2, j] = parse(Int, current_line[2])
    interface_info[3, j] = parse(Int, current_line[3])
    interface_info[4, j] = parse(Int, current_line[4])
    interface_info[5, j] = parse(Int, current_line[5])
    interface_info[6, j] = parse(Int, current_line[6])

    # count the number of physical boundaries
    if interface_info[4,j] == 0
      n_boundaries += 1
    end
    file_idx += 1
  end

  # work arrays to pull to correct corners of a given element (agnostic to curvature) and local
  # copies of the curved boundary information

  # readin an store the curved boundary information of the elements

  for j in 1:n_elements
    # pull the corner node IDs
    current_line           = split(file_lines[file_idx])
    element_node_ids[1, j] = parse(Int, current_line[1])
    element_node_ids[2, j] = parse(Int, current_line[2])
    element_node_ids[3, j] = parse(Int, current_line[3])
    element_node_ids[4, j] = parse(Int, current_line[4])
    for i in 1:4
      # pull the (x,y) values of these corners out of the nodes array
      cornerNodeVals[i, :] .= corner_nodes[:, element_node_ids[i, j]]
    end
    # pull the information to check if boundary is curved in order to read in additional data
    file_idx += 1
    current_line    = split(file_lines[file_idx])
    curved_check[1] = parse(Int, current_line[1])
    curved_check[2] = parse(Int, current_line[2])
    curved_check[3] = parse(Int, current_line[3])
    curved_check[4] = parse(Int, current_line[4])
    if sum(curved_check) == 0
      # quadrilateral element is straight sided
      element_is_curved[j] = false
      file_idx  += 1
      # read all the boundary names
      bndy_names[:, j] = split(file_lines[file_idx])
    else
      # quadrilateral element has at least one curved side
      element_is_curved[j] = true

      # flip node ordering to make sure the element is right-handed for the interpolations
      m1 = 1
      m2 = 2
      @views tempNodes[1, :] .= cornerNodeVals[4, :]
      @views tempNodes[2, :] .= cornerNodeVals[2, :]
      @views tempNodes[3, :] .= cornerNodeVals[3, :]
      @views tempNodes[4, :] .= cornerNodeVals[1, :]
      for i in 1:4
        if curved_check[i] == 0
          # when curved_check[i] is 0 then the "curve" from cornerNode(i) to cornerNode(i+1) is a
          # straight line. So we must construct the interpolant for this line
          for k in 1:mesh_nnodes
            curve_vals[k, 1] = tempNodes[m1, 1] + 0.5 * (cheby_nodes[k] + 1.0) * ( tempNodes[m2, 1]
                                                                                 - tempNodes[m1, 1])
            curve_vals[k, 2] = tempNodes[m1, 2] + 0.5 * (cheby_nodes[k] + 1.0) * ( tempNodes[m2, 2]
                                                                                 - tempNodes[m1, 2])
          end
        else
          # when curved_check[i] is 1 this curved boundary information is supplied by the mesh
          # generator. So we just read it into a work array
          for k in 1:mesh_nnodes
            file_idx      += 1
            current_line   = split(file_lines[file_idx])
            curve_vals[k, 1] = parse(RealT,current_line[1])
            curve_vals[k, 2] = parse(RealT,current_line[2])
          end
        end
        # construct the curve interpolant for the current side
        surface_curves[i, j] = CurvedSurfaceT(cheby_nodes, bary_weights, copy(curve_vals))
        # indexing update that contains a "flip" to ensure correct element orientation
        # if we need to construct the straight line "curves" when curved_check[i] == 0
        m1 += 1
        if i == 3
          m2 = 1
        else
          m2 += 1
        end
      end
      # finally read in the boundary names where "---" means an internal connection
      file_idx  += 1
      bndy_names[:, j] = split(file_lines[file_idx])
    end
    # one last increment to the global index to read the next piece of element information
    file_idx += 1
  end

  return n_boundaries
end

@inline Base.ndims(::UnstructuredQuadMesh) = 2
@inline Base.real(::UnstructuredQuadMesh{RealT}) where {RealT} = RealT

# Check if mesh is periodic
isperiodic(mesh::UnstructuredQuadMesh) = mesh.periodicity

Base.length(mesh::UnstructuredQuadMesh) = mesh.n_elements


function Base.show(io::IO, ::UnstructuredQuadMesh{RealT, CurvedSurfaceT}) where {RealT, CurvedSurfaceT}
  print(io, "UnstructuredQuadMesh{2, ", RealT, ", ", CurvedSurfaceT, "}")
end


function Base.show(io::IO, ::MIME"text/plain", mesh::UnstructuredQuadMesh{RealT, CurvedSurfaceT}) where {RealT, CurvedSurfaceT}
  if get(io, :compact, false)
    show(io, mesh)
  else
    summary_header(io, "UnstructuredQuadMesh{" * string(2) * ", " * string(RealT) * ", " * string(CurvedSurfaceT) * "}")
    summary_line(io, "mesh file", mesh.filename)
    summary_line(io, "number of elements", length(mesh))
    summary_line(io, "faces", mesh.n_surfaces)
    summary_line(io, "mesh polynomial degree", mesh.polydeg)
    summary_footer(io)
  end
end
