"""
    CurvedMesh{NDIMS, RealT<:Real} <: AbstractMesh{NDIMS}

A structured curved mesh.

Different numbers of cells per dimension are possible and arbitrary functions 
can be used as domain faces.

!!! warning "Experimental code"
    This mesh type is experimental and can change any time.
"""
mutable struct CurvedMesh{NDIMS, RealT<:Real} <: AbstractMesh{NDIMS}
  cells_per_dimension::NTuple{NDIMS, Int}
  mapping::Any # Not relevant for performance
  mapping_as_string::String
  periodicity::NTuple{NDIMS, Bool}
  current_filename::String
  unsaved_changes::Bool
end


"""
    CurvedMesh(cells_per_dimension, mapping, RealT; unsaved_changes=true, mapping_as_string=mapping2string(mapping, length(cells_per_dimension)))

Create a CurvedMesh of the given size and shape that uses `RealT` as coordinate type.

# Arguments
- `cells_per_dimension::NTupleE{NDIMS, Int}`: the number of cells in each dimension.
- `mapping`: a function of `NDIMS` variables to describe the mapping, which transforms 
             the reference mesh to the physical domain.
- `RealT::Type`: the type that should be used for coordinates.
- `periodicity`: either a `Bool` deciding if all of the boundaries are periodic or an `NTuple{NDIMS, Bool}` 
                 deciding for each dimension if the boundaries in this dimension are periodic.
- `unsaved_changes::Bool`: if set to `true`, the mesh will be saved to a mesh file.
- `mapping_as_string::String`: the code that defines the `mapping`.
                               If `CodeTracking` can't find the function definition, it can be passed directly here.
"""
function CurvedMesh(cells_per_dimension, mapping; RealT=Float64, periodicity=true, unsaved_changes=true, 
                    mapping_as_string=mapping2string(mapping, length(cells_per_dimension)))
  NDIMS = length(cells_per_dimension)

  # Convert periodicity to a Tuple of a Bool for every dimension
  if all(periodicity)
    # Also catches case where periodicity = true
    periodicity = ntuple(_->true, NDIMS)
  elseif !any(periodicity)
    # Also catches case where periodicity = false
    periodicity = ntuple(_->false, NDIMS)
  else
    # Default case if periodicity is an iterable
    periodicity = Tuple(periodicity)
  end

  return CurvedMesh{NDIMS, RealT}(Tuple(cells_per_dimension), mapping, mapping_as_string, periodicity, "", unsaved_changes)
end


"""
    CurvedMesh(cells_per_dimension, faces, RealT; unsaved_changes=true, faces_as_string=faces2string(faces))

Create a CurvedMesh of the given size and shape that uses `RealT` as coordinate type.

# Arguments
- `cells_per_dimension::NTupleE{NDIMS, Int}`: the number of cells in each dimension.
- `faces::NTuple{2*NDIMS}`: a tuple of `2 * NDIMS` functions that describe the faces of the domain.
                            Each function must take `NDIMS-1` arguments.
                            `faces[1]` describes the face onto which the face in negative x-direction 
                            of the unit hypercube is mapped. The face in positive x-direction of
                            the unit hypercube will be mapped onto the face described by `faces[2]`.
                            `faces[3:4]` describe the faces in positive and negative y-direction respectively 
                            (in 2D and 3D).
                            `faces[5:6]` describe the faces in positive and negative z-direction respectively (in 3D).
- `RealT::Type`: the type that should be used for coordinates.
- `periodicity`: either a `Bool` deciding if all of the boundaries are periodic or an `NTuple{NDIMS, Bool}` deciding for
                 each dimension if the boundaries in this dimension are periodic.
"""
function CurvedMesh(cells_per_dimension, faces::Tuple; RealT=Float64, periodicity=true)
  NDIMS = length(cells_per_dimension)

  validate_faces(faces)

  # Use the transfinite mapping with the correct number of arguments
  mapping = transfinite_mapping(faces)

  # Collect definitions of face functions in one string (separated by semicolons)
  face2substring(face) = code_string(face, ntuple(_ -> Float64, NDIMS-1))
  join_semicolon(strings) = join(strings, "; ")

  faces_definition = faces .|> face2substring .|> string |> join_semicolon

  # Include faces definition in `mapping_as_string` to allow for evaluation 
  # without knowing the face functions
  mapping_as_string = "$faces_definition; mapping = transfinite_mapping(faces)"

  return CurvedMesh(cells_per_dimension, mapping; RealT=RealT, periodicity=periodicity, mapping_as_string=mapping_as_string)
end


"""
    CurvedMesh(cells_per_dimension, coordinates_min, coordinates_max)

Create a CurvedMesh that represents a uncurved structured mesh with a rectangular domain.

# Arguments
- `cells_per_dimension::NTuple{NDIMS, Int}`: the number of cells in each dimension.
- `coordinates_min::NTuple{NDIMS, RealT}`: coordinate of the corner in the negative direction of each dimension.
- `coordinates_max::NTuple{NDIMS, RealT}`: coordinate of the corner in the positive direction of each dimension.
- `periodicity`: either a `Bool` deciding if all of the boundaries are periodic or an `NTuple{NDIMS, Bool}` deciding for
                 each dimension if the boundaries in this dimension are periodic.
"""
function CurvedMesh(cells_per_dimension, coordinates_min, coordinates_max; periodicity=true)
  NDIMS = length(cells_per_dimension)
  RealT = promote_type(eltype(coordinates_min), eltype(coordinates_max))

  mapping = coordinates2mapping(coordinates_min, coordinates_max)
  mapping_as_string = "coordinates_min = $coordinates_min; " *
                      "coordinates_max = $coordinates_max; " *
                      "mapping = coordinates2mapping(coordinates_min, coordinates_max)"
  return CurvedMesh(cells_per_dimension, mapping; RealT=RealT, periodicity=periodicity, mapping_as_string=mapping_as_string)
end


# Extract a string of the code that defines the mapping function
mapping2string(mapping, ndims) = string(code_string(mapping, ntuple(_ -> Float64, ndims)))


# Interpolate linearly between left and right value where s should be between -1 and 1
linear_interpolate(s, left_value, right_value) = 0.5 * ((1 - s) * left_value + (1 + s) * right_value)


# Convert min and max coordinates of a rectangle to the corresponding transformation mapping
function coordinates2mapping(coordinates_min::NTuple{1}, coordinates_max::NTuple{1})
  mapping(xi) = linear_interpolate(xi, coordinates_min[1], coordinates_max[1])
end

function coordinates2mapping(coordinates_min::NTuple{2}, coordinates_max::NTuple{2})
  mapping(xi, eta) = SVector(linear_interpolate(xi,  coordinates_min[1], coordinates_max[1]),
                             linear_interpolate(eta, coordinates_min[2], coordinates_max[2]))
end

function coordinates2mapping(coordinates_min::NTuple{3}, coordinates_max::NTuple{3})
  mapping(xi, eta, zeta) = SVector(linear_interpolate(xi,   coordinates_min[1], coordinates_max[1]),
                                   linear_interpolate(eta,  coordinates_min[2], coordinates_max[2]),
                                   linear_interpolate(zeta, coordinates_min[3], coordinates_max[3]))
end


# In 1D
# Linear mapping from the reference element to the domain described by the faces
function linear_mapping(x, faces)
  return linear_interpolate(x, faces[1](), faces[2]())
end


# In 2D
# Bilinear mapping from the reference element to the domain described by the faces
function bilinear_mapping(x, y, faces)
  x1 = faces[1](-1) # Bottom left
  x2 = faces[2](-1) # Bottom right
  x3 = faces[1](1) # Top left
  x4 = faces[2](1) # Top right

  return 0.25 * (x1 * (1 - x) * (1 - y) +
                 x2 * (1 + x) * (1 - y) +
                 x3 * (1 - x) * (1 + y) +
                 x4 * (1 + x) * (1 + y))
end


# In 3D
# Trilinear mapping from the reference element to the domain described by the faces
function trilinear_mapping(x, y, z, faces)
  x1 = faces[1](-1, -1) # mapped from (-1,-1,-1)
  x2 = faces[2](-1, -1) # mapped from ( 1,-1,-1)
  x3 = faces[1]( 1, -1) # mapped from (-1, 1,-1)
  x4 = faces[2]( 1, -1) # mapped from ( 1, 1,-1)
  x5 = faces[1](-1,  1) # mapped from (-1,-1, 1)
  x6 = faces[2](-1,  1) # mapped from ( 1,-1, 1)
  x7 = faces[1]( 1,  1) # mapped from (-1, 1, 1)
  x8 = faces[2]( 1,  1) # mapped from ( 1, 1, 1)

  return 0.125 * (x1 * (1 - x) * (1 - y) * (1 - z) +
                  x2 * (1 + x) * (1 - y) * (1 - z) +
                  x3 * (1 - x) * (1 + y) * (1 - z) +
                  x4 * (1 + x) * (1 + y) * (1 - z) +
                  x5 * (1 - x) * (1 - y) * (1 + z) +
                  x6 * (1 + x) * (1 - y) * (1 + z) +
                  x7 * (1 - x) * (1 + y) * (1 + z) +
                  x8 * (1 + x) * (1 + y) * (1 + z) )
end


# Use linear mapping in 1D
transfinite_mapping(faces::NTuple{2, Any}) = x -> linear_mapping(x, faces)

# In 2D
# Transfinite mapping from the reference element to the domain described by the faces
function transfinite_mapping(faces::NTuple{4, Any})
  mapping(x, y) = (linear_interpolate(x, faces[1](y), faces[2](y)) + 
                   linear_interpolate(y, faces[3](x), faces[4](x)) - 
                   bilinear_mapping(x, y, faces))
end


# In 3D
# Correction term for the Transfinite mapping
function correction_term_3d(x, y, z, faces)
  # Correction for x-terms
  c_x = linear_interpolate(x, linear_interpolate(y, faces[3](-1, z), faces[4](-1, z)) +
                              linear_interpolate(z, faces[5](-1, y), faces[6](-1, y)), 
                              linear_interpolate(y, faces[3]( 1, z), faces[4]( 1, z)) +
                              linear_interpolate(z, faces[5]( 1, y), faces[6]( 1, y)) )

  # Correction for y-terms
  c_y = linear_interpolate(y, linear_interpolate(x, faces[1](-1,  z), faces[2](-1,  z)) +
                              linear_interpolate(z, faces[5]( x, -1), faces[6]( x, -1)), 
                              linear_interpolate(x, faces[1]( 1,  z), faces[2]( 1,  z)) +
                              linear_interpolate(z, faces[5]( x,  1), faces[6]( x,  1)) )

  # Correction for x-terms
  c_z = linear_interpolate(z, linear_interpolate(x, faces[1](y, -1), faces[2](y, -1)) +
                              linear_interpolate(y, faces[3](x, -1), faces[4](x, -1)), 
                              linear_interpolate(x, faces[1](y,  1), faces[2](y,  1)) +
                              linear_interpolate(y, faces[3](x,  1), faces[4](x,  1)) )

  return 0.5 * (c_x + c_y + c_z)
end


# In 3D
# Transfinite mapping from the reference element to the domain described by the faces
function transfinite_mapping(faces::NTuple{6, Any})
  mapping(x, y, z) =  (linear_interpolate(x, faces[1](y, z), faces[2](y, z)) + 
                       linear_interpolate(y, faces[3](x, z), faces[4](x, z)) +
                       linear_interpolate(z, faces[5](x, y), faces[6](x, y)) -
                       correction_term_3d(x, y, z, faces) + 
                       trilinear_mapping(x, y, z, faces))
end


function validate_faces(faces::NTuple{2, Any}) end

function validate_faces(faces::NTuple{4, Any})
  @assert faces[1](-1) ≈ faces[3](-1) "faces[1](-1) needs to match faces[3](-1) (bottom left corner)"
  @assert faces[2](-1) ≈ faces[3]( 1) "faces[2](-1) needs to match faces[3](1) (bottom right corner)"
  @assert faces[1]( 1) ≈ faces[4](-1) "faces[1](1) needs to match faces[4](-1) (top left corner)"
  @assert faces[2]( 1) ≈ faces[4]( 1) "faces[2](1) needs to match faces[4](1) (top right corner)"
end

function validate_faces(faces::NTuple{6, Any})
  @assert (faces[1](-1, -1) ≈ 
           faces[3](-1, -1) ≈ 
           faces[5](-1, -1)) "faces[1](-1, -1), faces[3](-1, -1) and faces[5](-1, -1) need to match at (-1, -1, -1) corner"

  @assert (faces[2](-1, -1) ≈ 
           faces[3]( 1, -1) ≈ 
           faces[5]( 1, -1)) "faces[2](-1, -1), faces[3](1, -1) and faces[5](1, -1) need to match at (1, -1, -1) corner"
  
  @assert (faces[1]( 1, -1) ≈ 
           faces[4](-1, -1) ≈ 
           faces[5](-1,  1)) "faces[1](1, -1), faces[4](-1, -1) and faces[5](-1, 1) need to match at (-1, 1, -1) corner"
  
  @assert (faces[2]( 1, -1) ≈ 
           faces[4]( 1, -1) ≈ 
           faces[5]( 1,  1)) "faces[2](1, -1), faces[4](1, -1) and faces[5](1, 1) need to match at (1, 1, -1) corner"
  
  @assert (faces[1](-1,  1) ≈ 
           faces[3](-1,  1) ≈ 
           faces[6](-1, -1)) "faces[1](-1, 1), faces[3](-1, 1) and faces[6](-1, -1) need to match at (-1, -1, 1) corner"

  @assert (faces[2](-1,  1) ≈ 
           faces[3]( 1,  1) ≈ 
           faces[6]( 1, -1)) "faces[2](-1, 1), faces[3](1, 1) and faces[6](1, -1) need to match at (1, -1, 1) corner"
  
  @assert (faces[1]( 1,  1) ≈ 
           faces[4](-1,  1) ≈ 
           faces[6](-1,  1)) "faces[1](1, 1), faces[4](-1, 1) and faces[6](-1, 1) need to match at (-1, 1, 1) corner"
  
  @assert (faces[2]( 1,  1) ≈ 
           faces[4]( 1,  1) ≈ 
           faces[6]( 1,  1)) "faces[2](1, 1), faces[4](1, 1) and faces[6](1, 1) need to match at (1, 1, 1) corner"
end


# Check if mesh is periodic
isperiodic(mesh::CurvedMesh) = all(mesh.periodicity)
isperiodic(mesh::CurvedMesh, dimension) = mesh.periodicity[dimension]

@inline Base.ndims(::CurvedMesh{NDIMS}) where {NDIMS} = NDIMS
@inline Base.real(::CurvedMesh{NDIMS, RealT}) where {NDIMS, RealT} = RealT
Base.size(mesh::CurvedMesh) = mesh.cells_per_dimension
Base.size(mesh::CurvedMesh, i) = mesh.cells_per_dimension[i]
Base.axes(mesh::CurvedMesh) = map(Base.OneTo, mesh.cells_per_dimension)
Base.axes(mesh::CurvedMesh, i) = Base.OneTo(mesh.cells_per_dimension[i])


function Base.show(io::IO, ::CurvedMesh{NDIMS, RealT}) where {NDIMS, RealT}
  print(io, "CurvedMesh{", NDIMS, ", ", RealT, "}")
end


function Base.show(io::IO, ::MIME"text/plain", mesh::CurvedMesh{NDIMS, RealT}) where {NDIMS, RealT}
  if get(io, :compact, false)
    show(io, mesh)
  else
    summary_header(io, "CurvedMesh{" * string(NDIMS) * ", " * string(RealT) * "}")
    summary_line(io, "size", size(mesh))

    summary_line(io, "mapping", "")
    # Print code lines of mapping_as_string
    mapping_lines = split(mesh.mapping_as_string, ";")
    for i in eachindex(mapping_lines)
      summary_line(increment_indent(io), "line $i", strip(mapping_lines[i]))
    end
    summary_footer(io)
  end
end
