text:
# euler_funcs.jl
# this file contains all the functions that calculate values at a node (except
# for boundary related functions) as well as some additional helper functions.

#------------------------------------------------------------------------------
# function to calculate the Euler flux over the entire mesh
#------------------------------------------------------------------------------
@doc """
### EulerEquationMod.getEulerFlux

  This function calculates the Euler flux across the entire mesh by passing
  pieces of the eqn.q, eqn.aux_vars, eqn.f_xi and eqn.params to a low level
  function.  The flux is calculated in the xi and eta directions,
  scaled (mulitiplied) by the mapping jacobian (so that when performing the
  integral we don't have to explictly divide by the jacobian, it just cancels
  out with the jacobian factor introduced here.

  Calls writeFlux to do any requested output.

  This is a mid level function
"""->
# mid level function
function getEulerFlux(mesh::AbstractMesh{Tmsh},
                sbp::AbstractOperator,
                eqn::EulerData{Tsol, Tres, Tdim}, opts) where {Tmsh, Tsol, Tres, Tdim}
# calculate Euler flux in parametric coordinate directions, stores it in eqn.flux_parametric

  nrm = zeros(Tmsh, Tdim)
  for i=1:mesh.numEl  # loop over elements
    for j=1:mesh.numNodesPerElement  # loop over nodes on current element
      q_vals = ro_sview(eqn.q, :, j, i)
      aux_vars = ro_sview(eqn.aux_vars, :, j, i)
      # put this loop here (rather than outside) to we don't have to fetch
      # q_vals twice, even though writing to the flux vector is slower
      # it might be worth copying the normal vector rather than
      # doing an view

      for k=1:Tdim  # loop over dimensions
        for p=1:Tdim
          nrm[p] = mesh.dxidx[k, p, j, i]
        end
        # don't do an array view because strided views are type-unstable
#        nrm[1] = mesh.dxidx[k, 1, j, i]
#        nrm[2] = mesh.dxidx[k, 2, j, i]
        flux = sview(eqn.flux_parametric, :, j, i, k)

        # this will dispatch to the proper calcEulerFlux
        calcEulerFlux(eqn.params, q_vals, aux_vars, nrm, flux)
      end

    end
  end

  writeFlux(mesh, sbp, eqn, opts)

  return nothing
end

@doc """
### EulerEquationMod.writeFlux

  This function writes the real part of Euler flux to a file named Fxi.dat,
  space delimited, controlled by the input options 'writeflux', of type Bool.

  This is a high level function.
"""->
function writeFlux(mesh, sbp, eqn, opts)

   if !eqn.params.writeflux
     return nothing
   end

   fname = "Fxi.dat"
   rmfile(fname)
   writedlm(fname, real(eqn.flux_parametric))

   return nothing
end

@doc """
### EulerEquationMod.getEulerFlux2

  This function calcules the euler flux over the entire mesh directly (ie.
  does not call a low level function.  This function is deprecated, although
  useful for benchmarking purposes.  2D only.

  This is a mid level function
"""->
# this function is deprecated in favor of getEulerFlux()
# useful for benchmarking purposes
function getEulerFlux2(mesh::AbstractMesh{Tmsh}, sbp::AbstractOperator,
                       eqn::EulerData{Tsol}, opts) where {Tmsh, Tsol}
# calculates the Euler flux for every node in the xi and eta directions
# eqn is the equation type
# q is the 3D array (4 by nnodes per element by nel), of the conservative variables
# dxidx is the 4D array (2 by 2 x nnodes per element by nel) that specifies the direction of xi and eta in each element (output from mappingjacobian!)
# flux_parametric is populated with the flux in xi direction (same shape as q)
# F_eta is populated with flux in eta direction

# once the Julia developers fix slice notation and speed up subarrays, we won't have to
# vectorize like this (can calculate flux one node at a time inside a dedicated function

q = eqn.q
dxidx = mesh.dxidx
flux_parametric = sview(eqn.flux_parametric, :, :, :, 1)
F_eta = sview(eqn.flux_parametric, :, :, :, 2)

(ncomp, nnodes, nel) = size(q)  # get sizes of things

  for i=1:nel  # loop over elements
    for j=1:nnodes  # loop over nodes within element
      # get direction vector components (xi direction)
      nx = dxidx[1, 1, j, i]
      ny = dxidx[1, 2, j, i]
      # calculate pressure
      press = (eqn.params.gamma-1)*(q[4, j, i] - 0.5*(q[2, j, i]^2 + q[3, j, i]^2)/q[1, j, i])

      # calculate flux in xi direction
      # hopefully elements of q get stored in a register for reuse in eta direction
      U = (q[2, j, i]*nx + q[3, j, i]*ny)/q[1, j, i]
      flux_parametric[1, j, i] = q[1, j, i]*U
      flux_parametric[2, j, i] = q[2, j, i]*U + nx*press
      flux_parametric[3, j, i] = q[3, j, i]*U + ny*press
      flux_parametric[4, j, i] = (q[4, j, i] + press)*U

      # get direction vector components (eta direction)
      nx = dxidx[2, 1, j, i]
      ny = dxidx[2, 2, j, i]

      # calculate xi flux
      U = (q[2, j, i]*nx + q[3, j, i]*ny)/q[1, j, i]
      F_eta[1, j, i] = q[1, j, i]*U
      F_eta[2, j, i] = q[2, j, i]*U + nx*press
      F_eta[3, j, i] = q[3, j, i]*U + ny*press
      F_eta[4, j, i] = (q[4, j, i] + press)*U
    end
  end

  return nothing

end

"""
  Calculates the volume integrals for the weak form, computing the Euler flux
  as needed, rather than using eqn.flux_parametric

  Inputs:
    mesh
    sbp
    eqn: eqn.res is updated with the result
    opts
"""
function calcVolumeIntegrals_nopre(
           mesh::AbstractMesh{Tmsh},
           sbp::AbstractOperator,
           eqn::EulerData{Tsol, Tres, Tdim},
           opts) where {Tmsh, Tsol, Tres, Tdim}

  # flux in the parametric directions for a given element
  flux_el = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerElement, Tdim)
  nrm = zeros(Tmsh, mesh.dim)

  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      q_j = ro_sview(eqn.q, :, j, i)
      aux_vars_j = ro_sview(eqn.aux_vars, :, j, i)

      for k=1:Tdim
        flux_k = sview(flux_el, :, j, k)

        # get the direction vector
        for p=1:Tdim
          nrm[p] = mesh.dxidx[k, p, j, i]
        end
        # consider calculating all directions at once
        # not sure if that will help because the data dependencies are
        # really simple
        calcEulerFlux(eqn.params, q_j, aux_vars_j, nrm, flux_k)
      end  # end loop k
    end  # end loop j

    res_i = sview(eqn.res, :, :, i)
    for k=1:Tdim
      weakDifferentiateElement!(sbp, k, sview(flux_el, :, :, k), res_i, SummationByParts.Add(), true)
    end

  end  # end loop i

  return nothing
end  # end function

"""
  Calculate (S .*F)1 where S is the skew symmetric part of sbp.Q and F
  is a symmetric numerical flux function.  eqn.res is updated with the result.
  Methods are available for curvilinear and non-curvilinear meshes

  Inputs:
    mesh
    sbp
    eqn
    opts
    functor: the numerical flux function F, of type FluxType
"""
function calcVolumeIntegralsSplitForm(
                mesh::AbstractMesh{Tmsh},
                sbp::AbstractOperator,
                eqn::EulerData{Tsol, Tres, Tdim}, opts,
                functor::FluxType) where {Tmsh, Tsol, Tres, Tdim}

  if opts["use_staggered_grid"]
    # not planning on implementing the non-curvilinear version of this
    calcVolumeIntegralsSplitFormCurvilinear(mesh, mesh.mesh2, sbp, mesh.sbp2,
                                          eqn, opts, functor)
  else
    if mesh.coord_order == 1
      calcVolumeIntegralsSplitFormLinear(mesh, sbp, eqn, opts, functor)
    else
      calcVolumeIntegralsSplitFormCurvilinear(mesh, sbp, eqn, opts, functor)
    end
  end

  return nothing
end

@doc """
  Calculate (S .* F)1, where S is the skew-symmetric part of sbp.Q
  and F is a symmetric numerical flux function.  eqn.res is updated
  with the result.  Linear (non-curvilinear) meshes only
"""
function calcVolumeIntegralsSplitFormLinear(
                                        mesh::AbstractMesh{Tmsh},
                                        sbp::AbstractOperator,
                                        eqn::EulerData{Tsol, Tres, Tdim}, opts,
                                        functor::FluxType) where {Tmsh, Tsol, Tres, Tdim}

#  println("----- entered calcVolumeIntegralsSplitForm -----")
  dxidx = mesh.dxidx
  res = eqn.res
  q = eqn.q
  aux_vars = eqn.aux_vars
  params = eqn.params
  data = params.calc_volume_integrals_data
  @unpack data nrmD F_d S

  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      q_j = ro_sview(q, :, j, i)
      aux_vars_j = ro_sview(aux_vars, :, j, i)
      for k=1:(j-1)  # loop over lower triangle of S
        q_k = ro_sview(q, :, k, i)
        # calcaulate the normal vector in each parametric directions
        for d=1:Tdim
          # get the normal vector
          for p=1:Tdim
            nrmD[p, d] = dxidx[d, p, j, i]
          end
        end

        # calculate the numerical flux functions in all Tdim
        # directions at once
        functor(params, q_j, q_k, aux_vars_j, nrmD, F_d)

        @simd for d=1:Tdim
          # update residual
          @simd for p=1:(Tdim+2)
            res[p, j, i] -= 2*S[j, k, d]*F_d[p, d]
            res[p, k, i] += 2*S[j, k, d]*F_d[p, d]
          end

        end  # end d loop
      end  # end k loop

    end  # end j loop
  end  # end i loop

  return nothing
end

"""
  Calculate (S .* F)1, where S is the skew-symmetric part of sbp.Q
  and F is a symmetric numerical flux function.  eqn.res is updated
  with the result.  This function is used for curvilinear meshes.
"""
function calcVolumeIntegralsSplitFormCurvilinear(
                                        mesh::AbstractMesh{Tmsh},
                                        sbp::AbstractOperator,
                                        eqn::EulerData{Tsol, Tres, Tdim}, opts,
                                        functor::FluxType) where {Tmsh, Tsol, Tres, Tdim}

  dxidx = mesh.dxidx
  res = eqn.res
  q = eqn.q
  aux_vars = eqn.aux_vars
  params = eqn.params

  data = params.calc_volume_integrals_data
  @unpack data nrmD F_d Sx

  # S is calculated in x-y-z, so the normal vectors should be the unit normals
  fill!(nrmD, 0.0)
  for d=1:Tdim
    nrmD[d, d] = 1
  end

  for i=1:mesh.numEl
    # get S for this element
    dxidx_i = ro_sview(dxidx, :, :, :, i)
    calcSCurvilinear(sbp, dxidx_i, Sx)

    for j=1:mesh.numNodesPerElement
      q_j = ro_sview(q, :, j, i)
      aux_vars_j = ro_sview(aux_vars, :, j, i)
      for k=1:(j-1)  # loop over lower triangle of S
        q_k = ro_sview(q, :, k, i)

        # calculate the numerical flux functions in all Tdim
        # directions at once
        functor(params, q_j, q_k, aux_vars_j, nrmD, F_d)

        @simd for d=1:Tdim
          # update residual
          @simd for p=1:(Tdim+2)
            res[p, j, i] -= 2*Sx[j, k, d]*F_d[p, d]
            res[p, k, i] += 2*Sx[j, k, d]*F_d[p, d]
          end

        end  # end d loop
      end  # end k loop

    end  # end j loop

  end  # end i loop

  return nothing
end

"""
  This function does the interpolation from the solution grid to the flux
  grid for a single element

  **Inputs**:
   * params: a ParamType.  The qs_el1 and q_el1 fields are overwritten.
   * qs_el: the conservative variables for the element on the solution grid

  **Inputs/Outputs**:

   * aux_vars: array to be populated with the aux vars on the staggered grid
   * qf_el: conservative variables on the flux grid

"""
function interpolateElementStaggered(
                               params::ParamType{Tdim, :conservative},
                               mesh,
                               qs_el::AbstractMatrix,
                               aux_vars::AbstractMatrix,
                               qf_el::AbstractMatrix) where Tdim

  numNodesPerElement_s = size(qs_el, 2)
  numNodesPerElement_f = size(qf_el, 2)

  # temporary arrays
  @unpack params.interpolate_element_staggered_data wvars_s wvars_f

  for j=1:numNodesPerElement_s
    q_j = ro_sview(qs_el, :, j)
    w_j = sview(wvars_s, :, j)
    convertToIR(params, q_j, w_j)
  end

  # interpolate
  smallmatmat!(wvars_s, mesh.I_S2FT, wvars_f)

  # convert back to conservative
  for j=1:numNodesPerElement_f
    w_j = ro_sview(wvars_f, :, j)
    q_j = sview(qf_el, :, j)
    convertToConservativeFromIR_(params, w_j, q_j)
    aux_vars[1, j] = calcPressure(params, q_j)
  end

  return nothing
end

"""
  This function calculates I_S2F.'*(S .* F( uk(I_S2F wk), uk(I_S2F wk)))1.
  This is similar to the other method but for the staggered grid
  algorithm.

  Inputs:
    mesh_s: the mesh object for the solution grid
    mesh_f: the mesh object for the flux grid
    sbp_s: SBP operator for solution grid
    sbp_f: SBP operator for flux grid
    functor: the numerical flux functor ([`FluxType`](@ref)) used to compute F

  Inputs/Outputs:
    eqn: the equation object (implicitly on the solution grid).  eqn.res is
         updated with the result

  Aliasing Restrictions: none (the meshes and the sbp operators could alias,
                         in which case this algorithm reduces to the
                         non-staggered version
"""
function calcVolumeIntegralsSplitFormCurvilinear(
                                        mesh_s::AbstractMesh{Tmsh},
                                        mesh_f::AbstractMesh{Tmsh},
                                        sbp_s::AbstractOperator,
                                        sbp_f::AbstractOperator,
                                        eqn::EulerData{Tsol, Tres, Tdim}, opts,
                                        functor::FluxType) where {Tmsh, Tsol, Tres, Tdim}

  dxidx = mesh_f.dxidx
  res = eqn.res
  qf = eqn.q_flux
  params = eqn.params

  @unpack params.calc_volume_integrals_data nrmD F_d Sx res_f res_s
  aux_vars = zeros(Tres, 1, mesh_f.numNodesPerElement)

  # S is calculated in x-y-z, so the normal vectors should be the unit normals
  fill!(nrmD, 0.0)
  for d=1:Tdim
    nrmD[d, d] = 1
  end

  for i=1:mesh_f.numEl
    # get S for this element
    dxidx_i = ro_sview(dxidx, :, :, :, i)
    calcSCurvilinear(sbp_f, dxidx_i, Sx)

    fill!(res_f, 0.0)
    for j=1:mesh_f.numNodesPerElement
       q_j = ro_sview(qf, :, j, i)
       aux_vars[1, j] = calcPressure(eqn.params, q_j)
       aux_vars_j = ro_sview(aux_vars, :, j)
      for k=1:(j-1)  # loop over lower triangle of S
         q_k = ro_sview(qf, :, k, i)

        # calculate the numerical flux functions in all Tdim
        # directions at once
        functor(params, q_j, q_k, aux_vars_j, nrmD, F_d)

        @simd for d=1:Tdim
          # update residual
          @simd for p=1:(Tdim+2)
            res_f[p, j] -= 2*Sx[j, k, d]*F_d[p, d]
            res_f[p, k] += 2*Sx[j, k, d]*F_d[p, d]
          end

        end  # end d loop
      end  # end k loop

    end  # end j loop

    # reverse interpolate the residual
    #TODO: smallmatTmat might be faster
    # TODO: emulate BLAS interface to allow updating the rhs
    smallmatmat!(res_f, mesh_s.I_S2F, res_s)

    # accumulate into res
    @simd for j=1:mesh_s.numNodesPerElement
      @simd for k=1:mesh_s.numDofPerNode
        res[k, j, i] += res_s[k, j]
      end
    end

  end  # end i loop

  return nothing
end

# calculating the Euler flux at a node
#------------------------------------------------------------------------------
@doc """
### EulerEquationMod.calcEulerFlux

   This function calculates the Euler flux from the conservative variables at
   a single node in a particular direction.  2D only.

   Inputs:
   params  : ParamaterType{2, :conservative}
   q  : vector of conservative variables
   aux_vars : vector of auxiliary variables
   dir :  vector in direction to calculate the flux

   Inputs/Outputs:
   F  : vector to populate with the flux

   The Tdim paramater of params determine whether this method or the 3D
   version is called.

   This is a low level function
"""->
# low level function
function calcEulerFlux(params::ParamType{2, :conservative},
    q::AbstractArray{Tsol,1},
    aux_vars::AbstractArray{Tres, 1},
    dir::AbstractArray{Tmsh},  F::AbstractArray{Tsol,1}) where {Tmsh, Tsol, Tres}
# calculates the Euler flux in a particular direction at a point
# eqn is the equation type
# q is the vector (of length 4), of the conservative variables at the point
# aux_vars is the vector of auxiliary variables at the point
# dir is a vector of length 2 that specifies the direction
# F is populated with the flux (is a vector of length 4)
# 2D  only

  press = calcPressure(params, q)
#  press = getPressure(aux_vars)
#  press = @getPressure(aux_vars)
  U = (q[2]*dir[1] + q[3]*dir[2])/q[1]
  F[1] = q[1]*U
  F[2] = q[2]*U + dir[1]*press
  F[3] = q[3]*U + dir[2]*press
  F[4] = (q[4] + press)*U

  return nothing

end

@doc """
###EulerEquationMod.calcEulerFlux_revm

Compute the derivative of the euler flux in reverse mode w.r.t to unit vector
flux direction.

"""->
function calcEulerFlux_revm(params::ParamType{2, :conservative},
          q::AbstractArray{Tsol,1}, aux_vars,
          dir::AbstractArray{Tmsh,1},
          dir_bar::AbstractArray{Tmsh,1},
          F_bar::AbstractArray{Tres,1},
          ) where {Tmsh, Tsol, Tres}

  # Compute the reverse mode
  # Differentiate euler flux product with F_bar in reverse mode w.r.t dir to get
  # q_bar

  press = calcPressure(params, q)
  U = (q[2]*dir[1] + q[3]*dir[2])/q[1]

  # intermediate function that is only used in computing F so has to be reverse
  # diffed only in F_bar

  # dir_bar has dependence on both F and U
  # Reverse mode using F_bar
  dir_bar[1] += F_bar[2]*press
  dir_bar[2] += F_bar[3]*press
  U_bar = 0.0
  U_bar += F_bar[1]*q[1] + F_bar[2]*q[2] + F_bar[3]*q[3] + F_bar[4]*(q[4] + press)

  # Reverse mode using U_bar
  dir_bar[1] += U_bar*q[2]/q[1]
  dir_bar[2] += U_bar*q[3]/q[1]

  return nothing
end

function calcEulerFlux_revm(params::ParamType{3, :conservative},
            q::AbstractArray{Tsol,1}, aux_vars,
            dir::AbstractArray{Tmsh,1},
            dir_bar::AbstractArray{Tmsh,1},
            F_bar::AbstractArray{Tres,1}
            ) where {Tmsh,Tsol,Tres}

  # press = gami*(q[4] - 0.5*(q[2]^2 + q[3]^2 + q[4]^2)/q[1])
  press = calcPressure(params, q)

  U_bar = zero(Tres)
  # F[5] = (q[5] + press)*U
  U_bar += (q[5] + press)*F_bar[5]
  #F[4] = q[4]*U + dir[3]*press
  U_bar += q[4]*F_bar[4]
  dir_bar[3] += press*F_bar[4]
  # F[3] = q[3]*U + dir[2]*press
  U_bar += q[3]*F_bar[3]
  dir_bar[2] += press*F_bar[3]
  # F[2] = q[2]*U + dir[1]*press
  U_bar += q[2]*F_bar[2]
  dir_bar[1] += press*F_bar[2]
  # F[1] = q[1]*U
  U_bar += q[1]*F_bar[1]
  # U = (q[2]*dir[1] + q[3]*dir[2] + q[4]*dir[3])/q[1]
  dir_bar[1] += q[2]*U_bar/q[1]
  dir_bar[2] += q[3]*U_bar/q[1]
  dir_bar[3] += q[4]*U_bar/q[1]

  return nothing
end

@doc """
###EulerEquationMod.calcEulerFlux_revq

Compute the derivative of the Euler flux in reverse mode w.r.t q

"""->

function calcEulerFlux_revq(params::ParamType{2, :conservative},
          q::AbstractArray{Tsol,1}, q_bar::AbstractArray{Tsol, 1}, aux_vars,
          dir::AbstractArray{Tmsh,1}, F_bar::AbstractArray{Tres,1},
          ) where {Tmsh, Tsol, Tres}

  press = calcPressure(params, q)
  U = (q[2]*dir[1] + q[3]*dir[2])/q[1]

  U_bar = zero(Tsol)     # Initialize
  press_bar = zero(Tsol) #
  # Reverse diff F[4] = (q[4] + press)*U
  q_bar[4] += F_bar[4]*U
  U_bar += F_bar[4]*(q[4] + press)
  press_bar += F_bar[4]*U

  # Reverse diff F[3] = q[3]*U + dir[2]*press
  q_bar[3] += F_bar[3]*U
  U_bar += F_bar[3]*q[3]
  press_bar += F_bar[3]*dir[2]

  # Reverse diff F[2] =  q[2]*U + dir[1]*press
  q_bar[2] += F_bar[2]*U
  U_bar += F_bar[2]*q[2]
  press_bar += F_bar[2]*dir[1]

  # Reverse diff  F[1] = q[1]*U
  q_bar[1] += F_bar[1]*U
  U_bar += F_bar[1]*q[1]

  # Reverse diff U = (q[2]*dir[1] + q[3]*dir[2])/q[1]
  q_bar[2] += U_bar*dir[1]/q[1]
  q_bar[3] += U_bar*dir[2]/q[1]
  q_bar[1] -= U_bar*(q[2]*dir[1] + q[3]*dir[2])/(q[1]*q[1])

  # Reverse diff press = calcPressure(params, q)
  calcPressure_revq(params, q, q_bar, press_bar)

  return nothing
end

function calcEulerFlux_revq(params::ParamType{3, :conservative},
              q::AbstractArray{Tsol,1}, q_bar::AbstractArray{Tsol, 1}, aux_vars,
              dir::AbstractArray{Tmsh,1},
              F_bar::AbstractArray{Tsol,1},
              ) where {Tmsh, Tsol}

  #  Forward Sweep
  press = calcPressure(params, q)
  U = (q[2]*dir[1] + q[3]*dir[2] + q[4]*dir[3])/q[1]

  # Reverse sweep
  # F[5] = (q[5] + press)*U
  U_bar = F_bar[5]*(q[5] + press)
  press_bar = F_bar[5]*U
  q_bar[5] += F_bar[5]*U

  # F[4] = q[4]*U + dir[3]*press
  q_bar[4] += F_bar[4]*U
  U_bar += F_bar[4]*q[4]
  press_bar += F_bar[4]*dir[3]

  # F[3] = q[3]*U + dir[2]*press
  q_bar[3] += F_bar[3]*U
  U_bar += F_bar[3]*q[3]
  press_bar += F_bar[3]*dir[2]

  # F[2] = q[2]*U + dir[1]*press
  U_bar += F_bar[2]*q[2]
  q_bar[2] += F_bar[2]*U
  press_bar += F_bar[2]*dir[1]

  # F[1] = q[1]*U
  q_bar[1] += F_bar[1]*U
  U_bar += F_bar[1]*q[1]

  # U = (q[2]*dir[1] + q[3]*dir[2] + q[4]*dir[3])/q[1]
  q_bar[2] += U_bar*dir[1]/q[1]
  q_bar[3] += U_bar*dir[2]/q[1]
  q_bar[4] += U_bar*dir[3]/q[1]
  q_bar[1] -= U_bar*(q[2]*dir[1] + q[3]*dir[2] + q[4]*dir[3])/(q[1]*q[1])

  calcPressure_revq(params, q, q_bar, press_bar)

  return nothing
end

@doc """
# low level function
    Calculates the Euler flux from entropy variables

    Inputs:
    params : ParameterType{Tdim, :entropy}
    q : vector of entropy variables
    aux_vars : vector of auxiliary variables
    dir : vector specifying the direction to caculate the flux

    Inputs/Outputs:
    F  : vector to populate with the flux

    This is a low level function.  The static parameters of
    the ParameterType are used to dispatch to the right method for any
    combination of variable type or equation dimension.
"""->
function calcEulerFlux(params::ParamType{2, :entropy},
     q::AbstractArray{Tsol,1},
     aux_vars::AbstractArray{Tres, 1},
     dir::AbstractArray{Tmsh},  F::AbstractArray{Tsol,1}) where {Tmsh, Tsol, Tres}

  gamma = params.gamma
  gamma_1 = params.gamma_1

  # calculate some intermediate quantities
  k1 = 0.5*(q[2]^2 + q[3]^2)/q[4]  # a constant from Hughes' paper
  s = gamma - q[1] + k1    # entropy
    # internal energy (rho*i in Hughes) - not specific internal energy e
  rho_int = exp(-s/gamma_1)*(gamma_1/((-q[4])^gamma))^(1/gamma_1)
  U = q[2]*dir[1] + q[3]*dir[2]
  fac = rho_int/q[4]

  # now we can actually calculate the flux
  F[1] = q[4]*U*fac
  F[2] = (dir[1]*gamma_1*q[4] - q[2]*U)*fac
  F[3] = (dir[2]*gamma_1*q[4] - q[3]*U)*fac
  F[4] = U*(k1 - gamma)*fac

  return nothing
end

@doc """
### EulerEquationMod.calcEulerFlux
  This is the 3D method.  All arguments are same as the 2D version.
"""->
# low level function
function calcEulerFlux(params::ParamType{3, :conservative},
     q::AbstractArray{Tsol,1},
     aux_vars::AbstractArray{Tres, 1},
     dir::AbstractArray{Tmsh},
     F::AbstractArray{Tsol,1}) where {Tmsh, Tsol, Tres}
# calculates the Euler flux in a particular direction at a point
# eqn is the equation type
# q is the vector (of length 5), of the conservative variables at the point
# dir is a vector of length 3 that specifies the direction
# F is populated with the flux (is a vector of length 5)
# 3D  only

# once the Julia developers fix slice notation and speed up subarrays, we can make a faster
# vectorized version of this

  press = calcPressure(params, q)
  U = (q[2]*dir[1] + q[3]*dir[2] + q[4]*dir[3])/q[1]
  F[1] = q[1]*U
  F[2] = q[2]*U + dir[1]*press
  F[3] = q[3]*U + dir[2]*press
  F[4] = q[4]*U + dir[3]*press
  F[5] = (q[5] + press)*U

  return nothing

end

function calcEulerFlux(params::ParamType{3, :entropy},
     q::AbstractArray{Tsol,1},
     aux_vars::AbstractArray{Tres, 1},
     dir::AbstractArray{Tmsh},  F::AbstractArray{Tsol,1}) where {Tmsh, Tsol, Tres}

  gamma = params.gamma
  gamma_1 = params.gamma_1

  # calculate some intermediate quantities
  k1 = 0.5*(q[2]^2 + q[3]^2 + q[4]^2 )/q[5]  # a constant from Hughes' paper
  s = gamma - q[1] + k1    # entropy
    # internal energy (rho*i in Hughes) - not specific internal energy e
  rho_int = exp(-s/gamma_1)*(gamma_1/((-q[5])^gamma))^(1/gamma_1)
  U = q[2]*dir[1] + q[3]*dir[2] + q[4]*dir[3]
  fac = rho_int/q[5]

  # now we can actually calculate the flux
  F[1] = q[5]*U*fac
  F[2] = (dir[1]*gamma_1*q[5] - q[2]*U)*fac
  F[3] = (dir[2]*gamma_1*q[5] - q[3]*U)*fac
  F[4] = (dir[3]*gamma_1*q[5] - q[4]*U)*fac
  F[5] = U*(k1 - gamma)*fac

  return nothing
end

@doc """
### EulerEquationMod.getAuxVars

  This function calculates any extra variables that are stored across the mesh
  using the conservative variables eqn.q.  Currently only calculates pressure.

  Thi is a mid level function
"""->

#------------------------------------------------------------------------------
# functions for calculating additional quantities: pressure, entropy etc.
#------------------------------------------------------------------------------

# mid level function
function getAuxVars(mesh::AbstractMesh{Tmsh},
              eqn::EulerData{Tsol, Tres, Tdim}) where {Tmsh, Tsol, Tres, Tdim}
# calculate all auxiliary variables

  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      q_vals = ro_sview(eqn.q, :, j, i)

      # calculate pressure
      press = calcPressure(eqn.params, q_vals)
      @setPressure(eqn.aux_vars, j, i, press)
    end
  end

  return nothing
end

@doc """
### EulerEquationMod.calcPressure

  This function calculates the pressure from the conservative variables at a
  node in 2D.  It returns a single value.

  Inputs:
    params : ParamType{Tdim, var_type }
    q  : vector of conservative variables

  The parameter of params determines whether the 2D or 3D, conservative
  or entropy method is dispatched.

  This is a low level function.

  Aliasing restrictions: none
"""->
# low level function
function calcPressure(params::ParamType{2, :conservative},
                      q::AbstractArray{Tsol,1} ) where Tsol
  # calculate pressure for a node
  # q is a vector of length 4 of the conservative variables

  return  (params.gamma_1)*(q[4] - 0.5*(q[2]*q[2] + q[3]*q[3])/q[1])

end

@doc """
### EulerEquationMod.calcPressure
  This function calculates pressure using the entropy variables.

  Inputs:
    params : ParamType{2, :entropy}
    q  : vector of entropy varaibles

  returns pressure

  Aliasing restrictions: none
"""->
function calcPressure(params::ParamType{2, :entropy},
                      q::AbstractArray{Tsol,1}) where Tsol

  gamma = params.gamma
  gamma_1 = params.gamma_1

  k1 = 0.5*(q[2]*q[2] + q[3]*q[3])/q[4]  # a constant from Hughes' paper
  s = gamma - q[1] + k1    # entropy
  rho_int = exp(-s/gamma_1)*(gamma_1/((-q[4])^gamma))^(1/gamma_1)
  return gamma_1*rho_int
end

@doc """
### EulerEquationMod.calcPressure

  3D method.  See 2D method documentation
"""->
# low level function
function calcPressure(params::ParamType{3, :conservative},
                      q::AbstractArray{Tsol,1} ) where Tsol
  # calculate pressure for a node
  # q is a vector of length 5 of the conservative variables
  return  (params.gamma_1)*(q[5] - 0.5*(q[2]*q[2] + q[3]*q[3] + q[4]*q[4])/q[1])

end

function calcPressure(params::ParamType{3, :entropy},
                      q::AbstractArray{Tsol,1}) where Tsol

  gamma = params.gamma
  gamma_1 = params.gamma_1

  k1 = 0.5*(q[2]*q[2] + q[3]*q[3] + q[4]*q[4])/q[5]  # a constant from Hughes' paper
  s = gamma - q[1] + k1    # entropy
  rho_int = exp(-s/gamma_1)*(gamma_1/((-q[5])^gamma))^(1/gamma_1)
  return gamma_1*rho_int
end

@doc """
### EulerEquationMod.calcSpeedofSound

  This function calculates the speed of sound at a node and returns it.
  Methods are available for both conservative and entropy variables.

  Inputs:
    params:  ParamType{Tdim, var_type}
    q  vector of solution variables at a node

  Returns: speed of sound

  This is a low level function

  Aliasing restrictions: none
"""->

function calcSpeedofSound(params::ParamType{Tdim, :conservative},
                    q::AbstractArray{Tsol, 1}) where {Tdim, Tsol}
# calculates teh speed of sond at a node
  pressure = calcPressure(params, q)
  return sqrt((params.gamma*pressure)/q[1])

end

function calcSpeedofSound(params::ParamType{2, :entropy},
                          q::AbstractArray{Tsol, 1}) where Tsol
# calculate speed of sound using the same formula as conservative variables,
# just rewriting all variables in entropy variables

#  printbacktrace()
#  println("q = ", q)

  gamma = params.gamma
  gamma_1 = params.gamma_1
  k1 = 0.5*(q[2]^2 + q[3]^2)/q[4]  # a constant from Hughes' paper
  pressure = calcPressure(params, q)
  s = gamma - q[1] + k1    # entropy

  rho_int = exp(-s/gamma_1)*(gamma_1/((-q[4])^gamma))^(1/gamma_1)
  rho = -q[4]*rho_int

#  println("rho = ", rho)
  return sqrt((params.gamma*pressure)/rho)
end

function calcSpeedofSound(params::ParamType{3, :entropy},
                          q::AbstractArray{Tsol, 1}) where Tsol
# calculate speed of sound using the same formula as conservative variables,
# just rewriting all variables in entropy variables

#  printbacktrace()
#  println("q = ", q)

  gamma = params.gamma
  gamma_1 = params.gamma_1
  k1 = 0.5*(q[2]^2 + q[3]^2 + q[4]^2)/q[5]  # a constant from Hughes' paper
  pressure = calcPressure(params, q)
  s = gamma - q[1] + k1    # entropy

  rho_int = exp(-s/gamma_1)*(gamma_1/((-q[5])^gamma))^(1/gamma_1)
  rho = -q[5]*rho_int

#  println("rho = ", rho)
  return sqrt((params.gamma*pressure)/rho)
end

@doc """
### EulerEquationMod.calcEntropy

  This function calculates the entropy at a node and returns it.  Method are
  available for conservative and entropy variables, 2D or 3D

  Inputs:
    params: ParamType{Tdim, var_type}, used to dispatch to the right method.
    q: vector of solution variables at a node.

  Returns: entropy

  This is a low level function

  Aliasing Restrictions: none

"""->
function calcEntropy(params::ParamType{2, :conservative},
                     q::AbstractArray{Tsol,1} ) where Tsol

  gamma = params.gamma
  gamma_1 = params.gamma_1

  rho_int = q[4] - 0.5*(q[2]*q[2] + q[3]*q[3])/q[1]
  return log(gamma_1*rho_int/(q[1]^gamma))
end

function calcEntropy(params::ParamType{2, :entropy},
                     q::AbstractArray{Tsol,1}) where Tsol

  gamma = params.gamma
  gamma_1 = params.gamma_1

  return gamma - q[1] + 0.5*(q[2]*q[2] + q[3]*q[3])/q[4]
end

function calcEntropy(params::ParamType{3, :conservative},
                     q::AbstractArray{Tsol,1} ) where Tsol

  gamma = params.gamma
  gamma_1 = params.gamma_1

  rho_int = q[5] - 0.5*(q[2]*q[2] + q[3]*q[3] + q[4]*q[4])/q[1]
  return log(gamma_1*rho_int/(q[1]^gamma))
end

function calcEntropy(params::ParamType{3, :entropy},
                     q::AbstractArray{Tsol,1}) where Tsol

  gamma = params.gamma
  gamma_1 = params.gamma_1

  return gamma - q[1] + 0.5*(q[2]*q[2] + q[3]*q[3] + q[4]*q[4])/q[5]
end

"""
  This function calculates the entropy function U used to define the IR
  variablesat a node and returns it.   It does not return the physical entropy
  s.  This function is agnostic to the dimension of the equation.

  Inputs:
    params: a ParamType
    q: a vector of conservative variables at a node

  Outputs:
    U: the value of the entropy function

  Aliasing restrictions: none.
"""
function calcEntropyIR(params::ParamType{Tdim, :conservative},
               q::AbstractArray{Tsol,1} ) where {Tdim, Tsol}
# this calculate the entropy functiion U associated with the IR variables,
# not the physical entropy s

  gamma = params.gamma
  gamma_1 = params.gamma_1

  p = calcPressure(params, q)
  rho = q[1]
  U = -rho*(log(p) - gamma*log(rho))/gamma_1
  return U
end

"""
  This function calculates the vorticity at every node of an element
  (because vorticity requires derivatives of the velocity, it would be
  awkward to compute it at a node level).

  3D, conservative variables only

  Inputs:
    params: a ParamType
    q: numDofPerNode x numNodesPerElement array of conservative variables at
       each node
    dxidx: the 3 x 3 x numNodesPerElement  scaled  mapping jacobian at the node
    jac: numNodesPerElement vector of the mapping jacobian determinant at
         each node

  Input/Outputs:
    vorticity: a 3 x numNodesPerElement array  containing the vorticity in the
               x, y, and z directions at each node, overwritten

  Aliasing restrictions: from params: dxidx_element, velocities, velocity_deriv,
                                      velocity_deriv_xy
"""
function calcVorticity(params::ParamType{3, :conservative}, sbp,
     q::AbstractMatrix{Tsol}, dxidx::Abstract3DArray{Tmsh},
     jac::AbstractVector{Tmsh},
     vorticity::AbstractMatrix{Tres}) where {Tsol, Tmsh, Tres}

  Tdim = 3
  numDofPerNode = size(q, 1)
  numNodesPerElement = size(q, 2)
  data = params.calc_vorticity_data
  dxidx_unscaled = data.dxidx_unscaled
  velocities = data.velocities
  # velocity derivatives in parametric space
  # first dimension: velocity component, 3rd dimension parametric direction
  velocity_deriv = data.velocity_deriv
  fill!(velocity_deriv, 0.0)
  # velocity derivatives in xy space
  # first dimension: velocity, second dimension derivative direction
  velocity_deriv_xy = data.velocity_deriv_xy
  fill!(velocity_deriv_xy, 0.0)
  fill!(vorticity, 0.0)

  # unscale the mapping jacobian
  for i=1:numNodesPerElement
    jac_i = jac[i]
    for j=1:Tdim
      for k=1:Tdim
        dxidx_unscaled[k, j, i] = dxidx[k, j, i]*jac_i
      end
    end
  end

  # compute velocities
  for i=1:numNodesPerElement
    rho_i = q[1, i]
    for j=1:Tdim
      velocities[j, i] = q[j+1, i]/rho_i
    end
  end

  # differentiate velocities
  for d=1:Tdim
    differentiateElement!(sbp, d, velocities, sview(velocity_deriv, :, :, d))
  end

  for i=1:numNodesPerElement
    for v=1:3 # velocity component
      for cart_dim=1:3  # cartesian direction
        for para_dim=1:3  # parametric direction, summed
          velocity_deriv_xy[v, cart_dim, i] +=
              velocity_deriv[v, i, para_dim]*dxidx_unscaled[para_dim, cart_dim]
        end
      end
    end
  end

  # finally, compute the vorticity
  for i=1:numNodesPerElement
    vorticity[1, i] = velocity_deriv_xy[3, 2, i] - velocity_deriv_xy[2, 3, i]
    vorticity[2, i] = -velocity_deriv_xy[3, 1, i] + velocity_deriv_xy[1, 3, i]
    vorticity[3, i] = velocity_deriv_xy[2, 1, i] - velocity_deriv_xy[1, 2, i]
  end

  return nothing
end

#------------------------------------------------------------------------------
# function to calculate various coefficient matrices
#------------------------------------------------------------------------------

@doc """
### EulerEquationMod.calcA0

  This function calculates the A0 (ie. dq/dv, where q are the conservative
  and v are the entropy variables) for a node, and stores it A0

  The formation of A0 is given in Hughes

  Inputs:
    params : ParamType{Tdim, :entropy}
    q  : vector of entropy variables, length Tdim + 2

  Inputs/Outputs:
  A0 : (Tdim+2)x(Tdim+2) matrix populated with A0.  Overwritten

  Aliasing restrictions: none
"""->
function calcA0(params::ParamType{2, :entropy}, q::AbstractArray{Tsol,1},
                A0::AbstractArray{Tsol, 2}) where Tsol

  gamma = params.gamma
  gamma_1 = params.gamma_1

  k1 = 0.5*(q[2]^2 + q[3]^2)/q[4]  # a constant from Hughes' paper
  k2 = k1 - gamma
  k3 = k1*k1 - 2*gamma*k1 + gamma
#    k4 = k2 - gamma_1
  s = gamma - q[1] + k1    # entropy

  rho_int = exp(-s/gamma_1)*(gamma_1/((-q[4])^gamma))^(1/gamma_1)

  fac = rho_int/(gamma_1*q[4])

  # calculate the variables used in Hughes A.1
  c1 = gamma_1*q[4] - q[2]*q[2]
  c2 = gamma_1*q[4] - q[3]*q[3]

  d1 = -q[2]*q[3]

  e1 = q[2]*q[4]
  e2 = q[3]*q[4]

  # populate the matrix
  # the matrix is symmetric, but we don't use it because I think populating
  # the matrix will be faster if the matrix is write-only
  A0[1,1] = -q[4]*q[4]*fac
  A0[2,1] = e1*fac
  A0[3,1] = e2*fac
  A0[4,1] = q[4]*(1-k1)*fac
  A0[1,2] = e1*fac  # symmetric
  A0[2,2] = c1*fac
  A0[3,2] = d1*fac
  A0[4,2] = q[2]*k2*fac
  A0[1,3] = e2*fac  # symmetric
  A0[2,3] = d1*fac  # symmetric
  A0[3,3] = c2*fac
  A0[4,3] = q[3]*k2*fac
  A0[1,4] = q[4]*(1-k1)*fac  # symmetric
  A0[2,4] = q[2]*k2*fac   # symmetric
  A0[3,4] = q[3]*k2*fac  # symmetric
  A0[4,4] = -k3*fac

    return nothing
end

function calcA0(params::ParamType{3, :entropy}, q::AbstractArray{Tsol,1},
                A0::AbstractArray{Tsol, 2}) where Tsol

  gamma = params.gamma
  gamma_1 = params.gamma_1

  k1 = 0.5*(q[2]^2 + q[3]^2 + q[4]^2)/q[5]  # a constant from Hughes' paper
  k2 = k1 - gamma
  k3 = k1*k1 - 2*gamma*k1 + gamma
  s = gamma - q[1] + k1    # entropy

  rho_int = exp(-s/gamma_1)*(gamma_1/((-q[5])^gamma))^(1/gamma_1)
  fac = rho_int/(gamma_1*q[5])

  # calculate the variables used in Hughes A.1
  c1 = gamma_1*q[5] - q[2]*q[2]
  c2 = gamma_1*q[5] - q[3]*q[3]
  c3 = gamma_1*q[5] - q[4]*q[4]
  d1 = -q[2]*q[3]
  d2 = -q[2]*q[4]
  d3 = -q[3]*q[4]

  e1 = q[2]*q[5]
  e2 = q[3]*q[5]
  e3 = q[4]*q[5]

  # populate the matrix
  # the matrix is symmetric, but we don't use it because I think populating
  # the matrix will be faster if the matrix is write-only
  A0[1,1] = -q[5]*q[5]*fac
  A0[2,1] = e1*fac
  A0[3,1] = e2*fac
  A0[4,1] = e3*fac
  A0[5,1] = q[5]*(1-k1)*fac
  A0[1,2] = e1*fac  # symmetric
  A0[2,2] = c1*fac
  A0[3,2] = d1*fac
  A0[4,2] = d2*fac
  A0[5,2] = q[2]*k2*fac
  A0[1,3] = e2*fac  # symmetric
  A0[2,3] = d1*fac  # symmetric
  A0[3,3] = c2*fac
  A0[4,3] = d3*fac
  A0[5,3] = q[3]*k2*fac
  A0[1,4] = e3*fac  # symmetric
  A0[2,4] = d2*fac  # symmetric
  A0[3,4] = d3*fac  # symmetric
  A0[4,4] = c3*fac
  A0[5,4] = q[4]*k2*fac
  A0[1,5] = q[5]*(1-k1)*fac  # symmetric
  A0[2,5] = q[2]*k2*fac  # symmetric
  A0[3,5] = q[3]*k2*fac  # symmetric
  A0[4,5] = q[4]*k2*fac  # symmetric
  A0[5,5] = -k3*fac

  return nothing
end

@doc """
# EulerEquationMod.calcA0Inv

  Calculates inv(A0), where A0 = dq/dv, where q are the conservative variables
  at a node and v are the entropy varaibles at a node, using the entropy
  variables.

  Inputs:
    params : ParamType{Tdim, :entropy}
    q  : vector (length 4 or 5) of entropy variables at a node

  Inputs/Outputs:
    A0inv : matrix to be populated with inv(A0).  Overwritten.

  Aliasing restrictions: none
"""->
function calcA0Inv(params::ParamType{2, :entropy},
             q::AbstractArray{Tsol,1},
             A0inv::AbstractArray{Tsol, 2}) where Tsol
  gamma = params.gamma
  gamma_1 = params.gamma_1

  k1 = 0.5*(q[2]^2 + q[3]^2)/q[4]  # a constant from Hughes' paper
   s = gamma - q[1] + k1    # entropy

  rho_int = exp(-s/gamma_1)*(gamma_1/((-q[4])^gamma))^(1/gamma_1)

  fac = -1/(rho_int*q[4])

  d1 = -q[2]*q[3]
  e1 = q[2]*q[4]
  e2 = q[3]*q[4]

  A0inv[1,1] = (k1*k1 + gamma)*fac
  A0inv[2,1] = k1*q[2]*fac
  A0inv[3,1] = k1*q[3]*fac
  A0inv[4,1] = (k1 + 1)*q[4]*fac
  A0inv[1,2] = k1*q[2]*fac  # symmetry
  A0inv[2,2] = (q[2]*q[2] - q[4])*fac
  A0inv[3,2] = -d1*fac
  A0inv[4,2] = e1*fac
  A0inv[1,3] = k1*q[3]*fac  # symmetry
  A0inv[2,3] = -d1*fac  # symmetry
  A0inv[3,3] = (q[3]*q[3] - q[4])*fac
  A0inv[4,3] = e2*fac
  A0inv[1,4] = (k1 + 1)*q[4]*fac  # symmetry
  A0inv[2,4] = e1*fac  # symmetry
  A0inv[3,4] = e2*fac  # symmetry
  A0inv[4,4] = q[4]*q[4]*fac

    return nothing
end

function calcA0Inv(params::ParamType{3, :entropy},
             q::AbstractArray{Tsol,1},
             A0inv::AbstractArray{Tsol, 2}) where Tsol
  gamma = params.gamma
  gamma_1 = params.gamma_1

  k1 = 0.5*(q[2]^2 + q[3]^2 + q[4]^2)/q[5]  # a constant from Hughes' paper
   s = gamma - q[1] + k1    # entropy

  rho_int = exp(-s/gamma_1)*(gamma_1/((-q[5])^gamma))^(1/gamma_1)

  fac = -1/(rho_int*q[5])

  d1 = -q[2]*q[3]
  d2 = -q[2]*q[4]
  d3 = -q[3]*q[4]

  e1 = q[2]*q[5]
  e2 = q[3]*q[5]
  e3 = q[4]*q[5]

  A0inv[1,1] = (k1*k1 + gamma)*fac
  A0inv[2,1] = k1*q[2]*fac
  A0inv[3,1] = k1*q[3]*fac
  A0inv[4,1] = k1*q[4]*fac
  A0inv[5,1] = (k1 + 1)*q[5]*fac
  A0inv[1,2] = k1*q[2]*fac  # symmetry
  A0inv[2,2] = (q[2]*q[2] - q[5])*fac
  A0inv[3,2] = -d1*fac
  A0inv[4,2] = -d2*fac
  A0inv[5,2] = e1*fac
  A0inv[1,3] = k1*q[3]*fac  # symmetry
  A0inv[2,3] = -d1*fac  # symmetry
  A0inv[3,3] = (q[3]*q[3] - q[5])*fac
  A0inv[4,3] = -d3*fac
  A0inv[5,3] = e2*fac
  A0inv[1,4] = k1*q[4]*fac  # symmetry
  A0inv[2,4] = -d2*fac  # symmetry
  A0inv[3,4] = -d3*fac  # symmetry
  A0inv[4,4] = (q[4]*q[4] - q[5])*fac
  A0inv[5,4] = e3*fac
  A0inv[1,5] = (k1 + 1)*q[5]*fac  # symmetry
  A0inv[2,5] = e1*fac  # symmetry
  A0inv[3,5] = e2*fac  # symmetry
  A0inv[4,5] = e3*fac  # symmetry
  A0inv[5,5] = q[5]*q[5]*fac

    return nothing
end

@doc """
### EulerEquationMod.calcA0

  This function calculates the A0 (ie. dq/dq, where q are the conservative
  variables at a node), and stores it A0.  This function is provided for
  compatability purposes

  Inputs:
    params : ParamType{2, :entropy}
    q  : vector of conservative variables, length 4

  Inputs/Outputs:
  A0 : 4x4 (or 5x5 in 3D)  matrix populated with A0.  Overwritten

  Aliasing restrictions: none
"""->
function calcA0(params::ParamType{Tdim, :conservative}, q::AbstractArray{Tsol,1},
          A0::AbstractArray{Tsol, 2}) where {Tdim, Tsol}

  for i=1:length(q)
    A0[i,i] = one(Tsol)
  end

  return nothing
end

@doc """
# EulerEquationMod.calcA0Inv

  Calculates inv(A0), where A0 = dq/dq, where q are the conservative variables
  at a node.  This function is provided for compatability purposes

  Inputs:
    params : ParamType{2, :entropy}
    q  : vector (length Tdim + 2) of conservative variables at a node

  Inputs/Outputs:
    A0inv : matrix to be populated with inv(A0).  Overwritten.

  Aliasing restrictions: none
"""->
function calcA0Inv(params::ParamType{Tdim, :conservative},
       q::AbstractArray{Tsol,1},
       A0inv::AbstractArray{Tsol, 2}) where {Tdim, Tsol}

  for i=1:length(q)
    A0inv[i,i] = one(Tsol)
  end

  return nothing
end

@doc """
### EulerEquationMod.matVecA0inv

  This function takes a 3D array and multiplies it in place by the inv(A0)
  matrix (calculated at each node), inplace, (where A0 =dq/dv, where q are the
  conservative variables and v are some other variables), inplace.
  Methods are available for conservative and entropy variables.

  For conservative variables, A0 is the identity matrix and this is a no-op

  Inputs:
    mesh
    sbp
    eqn
    opts

  Inputs/Outputs:
    res_arr: the array to multiply

  Aliasing restrictions: none
"""->
function matVecA0inv(mesh::AbstractMesh{Tmsh},
                     sbp::AbstractOperator,
                     eqn::EulerData{Tsol, Tres, Tdim, :entropy}, opts,
                     res_arr::AbstractArray{Tsol, 3}) where {Tmsh, Tsol, Tdim, Tres}
# multiply a 3D array by inv(A0) in-place, useful for explicit time stepping
# res_arr *can* alias eqn.q safely
  A0inv = Array{Tsol}(mesh.numDofPerNode, mesh.numDofPerNode)
  res_vals = Array{Tsol}(mesh.numDofPerNode)
  q_vals = Array{Tsol}(mesh.numDofPerNode)
  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      # copy values into workvec
      for k=1:mesh.numDofPerNode
        q_vals[k] = eqn.q[k, j, i]
        res_vals[k] = res_arr[k, j, i]
      end

      res_view = sview(res_arr, :, j, i)
      # get A0Inv for this node
      calcA0Inv(eqn.params, q_vals, A0inv)

      smallmatvec!(A0inv, res_vals, res_view)
    end
  end

  return nothing
end

#TODO: move this to Utils?
# no-op, because for conservative variables this is A0inv is the identity matrix
function matVecA0inv(mesh::AbstractMesh{Tmsh},
                     sbp::AbstractOperator,
                     eqn::AbstractSolutionData,
                     opts, res_arr::AbstractArray{Tsol, 3}) where {Tmsh, Tsol}

  return nothing
end

@doc """
### EulerEquationMod.matVecA0

  This function is the same as matVecA0inv, except it multiplies by A0 not
  A0inv.  See its documention.

"""->
function matVecA0(mesh::AbstractMesh{Tmsh},
sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim, :entropy},
opts, res_arr::AbstractArray{Tsol, 3}) where {Tmsh, Tsol, Tdim, Tres}
# multiply a 3D array by inv(A0) in-place, useful for explicit time stepping
# res_arr *can* alias eqn.q safely
# a non-alias tolerant implemention wold avoid copying q_vals
  A0 = Array{Tsol}(mesh.numDofPerNode, mesh.numDofPerNode)
  res_vals = Array{Tsol}(mesh.numDofPerNode)
  q_vals = Array{Tsol}(mesh.numDofPerNode)
  for i=1:mesh.numEl
    for j=1:mesh.numNodesPerElement
      # copy values into workvec
      for k=1:mesh.numDofPerNode
        q_vals[k] = eqn.q[k, j, i]
        res_vals[k] = res_arr[k, j, i]
      end

      res_view = sview(res_arr, :, j, i)
      # get A0Inv for this node
      calcA0(eqn.params, q_vals, A0)

      smallmatvec!(A0, res_vals, res_view)
    end
  end

  return nothing
end

#no-op
function matVecA0(mesh::AbstractMesh{Tmsh}, sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim, :conservative}, opts, res_arr::AbstractArray{Tsol, 3}) where {Tmsh, Tsol, Tdim, Tres}

  return nothing
end

@doc """
### EulerEquationMod.calcA1

  This function calculates the A1 (ie. dF1/dq, where F1 is the first column of
  the Euler flux) for a node, aka the flux
  Jacobian of the Euler flux in the x direction.  Methods are available for
  both conservative and entropy variables.

  The formation of A1 is given in Hughes paper

  Inputs:
    params : ParamType{2, :entropy}
    q  : vector of entropy variables, length 4

  Inputs/Outputs:
  A1 : 4x4 matrix to be populated.  Overwritten

"""->

function calcA1(params::ParamType{2, :conservative},
                q::AbstractArray{Tsol,1}, A1::AbstractArray{Tsol, 2}) where Tsol
  gamma_1 = params.gamma_1
  gamma = params.gamma
  u = q[2]/q[1] # Get velocity in the x-direction
  v = q[3]/q[1] # Get velocity in the x-direction

  intvar = gamma_1*(u*u + v*v)/2
  a1 = intvar*2 - gamma*q[4]/q[1]
  A1[1,1] = 0
  A1[2,1] = intvar - u*u
  A1[3,1] = -u*v
  A1[4,1] = a1*u

  A1[1,2] = 1
  A1[2,2] = (3 - gamma)*u
  A1[3,2] = v
  A1[4,2] = gamma*q[4]/q[1] - intvar - gamma_1*u*u

  A1[1,3] = 0
  A1[2,3] = -gamma_1*v
  A1[3,3] = u
  A1[4,3] = -gamma_1*u*v

  A1[1,4] = 0
  A1[2,4] = gamma_1
  A1[3,4] = 0
  A1[4,4] = gamma*u

  return nothing
end

function calcA1(params::ParamType{2, :entropy}, q::AbstractArray{Tsol,1},
                A1::AbstractArray{Tsol, 2}) where Tsol

  gamma = params.gamma
  gamma_1 = params.gamma_1

  k1 = 0.5*(q[2]^2 + q[3]^2)/q[4]  # a constant from Hughes' paper
  k2 = k1 - gamma
  k3 = k1*k1 - 2*gamma*k1 + gamma
  k4 = k2 - gamma_1
  k5 = k2*k2 - gamma_1*(k1 + k2)

  s = gamma - q[1] + k1    # entropy

  rho_int = exp(-s/gamma_1)*(gamma_1/((-q[4])^gamma))^(1/gamma_1)

  fac = rho_int/(gamma_1*q[4]*q[4])

  # calculate the variables used in Hughes A.1
  c1 = gamma_1*q[4] - q[2]*q[2]
  c2 = gamma_1*q[4] - q[3]*q[3]

  d1 = -q[2]*q[3]

  e1 = q[2]*q[4]
#  e2 = q[3]*q[4]

  # populate the matrix
  # the matrix is symmetric, but we don't use it because I think populating
  # the matrix will be faster if the matrix is write-only
  A1[1,1] = e1*q[4]*fac
  A1[2,1] = c1*q[4]*fac
  A1[3,1] = d1*q[4]*fac
  A1[4,1] = k2*e1*fac
  A1[1,2] = c1*q[4]*fac  # symmetric
  A1[2,2] = -(c1 + 2*gamma_1*q[4])*q[2]*fac
  A1[3,2] = -c1*q[3]*fac
  A1[4,2] = (c1*k2 + gamma_1*q[2]*q[2])*fac
  A1[1,3] = d1*q[4]*fac  # symmetric
  A1[2,3] = -c1*q[3]*fac  # symmetric
  A1[3,3] = -c2*q[2]*fac
  A1[4,3] = k4*d1*fac
  A1[1,4] = k2*e1*fac  # symmetric
  A1[2,4] = (c1*k2 + gamma_1*q[2]*q[2])*fac   # symmetric
  A1[3,4] = k4*d1*fac  # symmetric
  A1[4,4] = k5*q[2]*fac

    return nothing
end

@doc """
### EulerEquationMod.calcA2

  This function calculates A2 (ie. dF2/dq, where F2 is the second column of the
  Euler flux, aka the flux jacobian in the y direction.
  Methods are available for both conservative and entropy variables.

  The formation of A2 is given in Hughes

  Inputs:
    params : ParamType{2, :entropy},
    q  : vector of entropy variables, length 4
  Inputs/Outputs:
  A2 : 4x4 matrix to be populated.  Overwritten

"""->
function calcA2(params::ParamType{2, :conservative},
                q::AbstractArray{Tsol,1}, A2::AbstractArray{Tsol, 2}) where Tsol
  gamma_1 = params.gamma_1
  gamma = params.gamma
  u = q[2]/q[1] # Get velocity in the x-direction
  v = q[3]/q[1] # Get velocity in the x-direction

  intvar = gamma_1*(u*u + v*v)/2
  a1 = intvar*2 - gamma*q[4]/q[1]

  A2[1,1] = 0
  A2[2,1] = -u*v
  A2[3,1] = intvar - v*v
  A2[4,1] = a1*v
  A2[1,2] = 0
  A2[2,2] = v
  A2[3,2] = -gamma_1*u
  A2[4,2] = -gamma_1*u*v
  A2[1,3] = 1
  A2[2,3] = u
  A2[3,3] = (3 - gamma)*v
  A2[4,3] = gamma*q[4]/q[1] - intvar - gamma_1*v*v
  A2[1,4] = 0
  A2[2,4] = 0
  A2[3,4] = gamma_1
  A2[4,4] = gamma*v

  return nothing
end

function calcA2(params::ParamType{2, :entropy}, q::AbstractArray{Tsol,1},
                A2::AbstractArray{Tsol, 2}) where Tsol

  gamma = params.gamma
  gamma_1 = params.gamma_1

  k1 = 0.5*(q[2]^2 + q[3]^2)/q[4]  # a constant from Hughes' paper
  k2 = k1 - gamma
  k3 = k1*k1 - 2*gamma*k1 + gamma
  k4 = k2 - gamma_1
  k5 = k2*k2 - gamma_1*(k1 + k2)

  s = gamma - q[1] + k1    # entropy

  rho_int = exp(-s/gamma_1)*(gamma_1/((-q[4])^gamma))^(1/gamma_1)

  fac = rho_int/(gamma_1*q[4]*q[4])

  # calculate the variables used in Hughes A.1
  c1 = gamma_1*q[4] - q[2]*q[2]
  c2 = gamma_1*q[4] - q[3]*q[3]

  d1 = -q[2]*q[3]

#  e1 = q[2]*q[4]
  e2 = q[3]*q[4]

  # populate the matrix
  # the matrix is symmetric, but we don't use it because I think populating
  # the matrix will be faster if the matrix is write-only
  A2[1,1] = e2*q[4]*fac
  A2[2,1] = d1*q[4]*fac
  A2[3,1] = c2*q[4]*fac
  A2[4,1] = k2*e2*fac
  A2[1,2] = d1*q[4]*fac  # symmetric
  A2[2,2] = -c1*q[3]*fac
  A2[3,2] = -c2*q[2]*fac
  A2[4,2] = k4*d1*fac
  A2[1,3] = c2*q[4]*fac  # symmetric
  A2[2,3] = -c2*q[2]*fac  # symmetric
  A2[3,3] = -(c2 + 2*gamma_1*q[4])*q[3]*fac
  A2[4,3] = (c2*k2 + gamma_1*q[3]*q[3])*fac
  A2[1,4] = k2*e2*fac  # symmetric
  A2[2,4] = k4*d1*fac   # symmetric
  A2[3,4] = (c2*k2 + gamma_1*q[3]*q[3])*fac  # symmetric
  A2[4,4] = k5*q[3]*fac

    return nothing
end

@doc """

"""->

function calcSteadyFluxJacobian(params, A, q, nrm)

  u = q[2]/q[1]
  v = q[3]/q[1]

  a1 = params.gamma*q[4]/q[1]
  theta = nrm[1]*u + nrm[2]*v
  phi_2 = 0.5*params.gamma_1*(u*u + v*v)

  A[1,1] = 0.0
  A[2,1] = -u*theta + nrm[1]*phi_2
  A[3,1] = -v*theta + nrm[2]*phi_2
  A[4,1] = theta*(phi_2-a1)

  A[1,2] = nrm[1]
  A[2,2] = theta - (params.gamma - 2)*nrm[1]*u
  A[3,2] = nrm[1]*v - params.gamma_1*nrm[2]*u
  A[4,2] = nrm[1]*a1 - params.gamma_1*u*theta

  A[1,3] = nrm[2]
  A[2,3] = nrm[2]*u - params.gamma_1*nrm[1]*v
  A[3,3] = theta - (params.gamma - 2)*nrm[2]*v
  A[4,3] = nrm[2]*a1 - params.gamma_1*v*theta

  A[1,4] = 0.0
  A[2,4] = params.gamma_1*nrm[1]
  A[3,4] = params.gamma_1*nrm[2]
  A[4,4] = params.gamma*theta

  return nothing
end

@doc """
### EulerEquationMod.calcMaxWaveSpeed

  This function calculates the maximum wave speed (ie. acoustic wave speed)
  present in the domain and returns it.  Methods are available for conservative   and entropy variables.  This function uses eqn.q_vec, not eqn.q to do the
  calculation.

  This is a mid level function
"""->
function calcMaxWaveSpeed(mesh, sbp,
        eqn::EulerData{Tsol, Tres, Tdim, :conservative}, opts) where {Tsol, Tdim, Tres}
# calculate the maximum wave speed (ie. characteristic speed) on the mesh
# uses solution vector q, not array
  q = eqn.q_vec
  max_speed = zero(Float64)
  for i=1:mesh.numDofPerNode:length(q)
    q_i = ro_sview(q, i:(i+mesh.numDofPerNode - 1))
    a = calcSpeedofSound(eqn.params, q_i)
    u_nrm = zero(Tsol)
    for j=1:Tdim
      u_j = q_i[j+1]/q_i[1]
      u_nrm += u_j*u_j
    end
    u_norm = sqrt(real(u_nrm))
    speed = real(a) + u_norm

    if speed > max_speed
      max_speed = speed
    end  # end if statement
  end  # end loop over vector q

  max_speed = MPI.Allreduce(max_speed, MPI.MAX, mesh.comm)
  return max_speed
end  # end function

function calcMaxWaveSpeed(mesh, sbp,
              eqn::EulerData{Tsol, Tres, 2, :entropy}, opts) where {Tsol, Tres}
# calculate the maximum wave speed (ie. characteristic speed) on the mesh
# uses solution vector q, not array
  q = eqn.q_vec
  gamma = eqn.params.gamma
  gamma_1 = eqn.params.gamma_1

  max_speed = zero(eltype(q))
  for i=1:mesh.numDofPerNode:length(q)
    q_i = ro_sview(q, i:(i+mesh.numDofPerNode - 1))
    a = calcSpeedofSound(eqn.params, q_i)

    # this is dimension specific
    k1 = 0.5*(q_i[2]*q_i[2] + q_i[3]*q_i[3])/q_i[4]  # a constant from Hughes' paper
    s = gamma - q_i[1] + k1    # entropy
    rho_int = exp(-s/gamma_1)*(gamma_1/((-q_i[4])^gamma))^(1/gamma_1)
    rho = -rho_int*q_i[4]
    u_norm = zero(Tsol)
    for j=1:2
      u_j = rho_int*q[j+1]/rho
      u_norm += u_j*u_j
    end
    u_norm = sqrt(u_norm)
    speed = a + u_norm

    if speed > max_speed
      max_speed = speed
    end  # end if statement
  end  # end loop over vector q

  return max_speed
end  # end function

"""
  Calculate the maximum wave speed for a given state

  **Inputs**

   * params: a ParamType
   * qL: a vector of conservative variables at a node
   * dir: a direction vector, length 2 in 2D and 3 in 3D

  **Outputs**

   * lambda_max: the maximum wave speed

  Aliasing restrictions: none
"""
function getLambdaMax(params::ParamType{Tdim},
                      qL::AbstractVector{Tsol},
                      dir::AbstractVector{Tmsh}) where {Tsol, Tmsh, Tdim}

  Tres = promote_type(Tsol, Tmsh)
  gamma = params.gamma
  Un = zero(Tres)
  dA = zero(Tmsh)
  rhoLinv = 1/qL[1]

  pL = calcPressure(params, qL)
  aL = sqrt(gamma*pL*rhoLinv)  # speed of sound

  for i=1:Tdim
    Un += dir[i]*qL[i+1]*rhoLinv
    dA += dir[i]*dir[i]
  end

  dA = sqrt(dA)

  lambda_max = absvalue3(Un) + dA*aL

  return lambda_max
end

"""
  Method to compute the max wave speed in the volume (ie not in any
  particular direction: sqrt(u*u + w*w) + sqrt(gamma*p/rho)

  **Inputs**

   * params: ParamType
   * qL; vector of conservative variables at a node

  **Outputs**

   * lambda_max
"""
function getLambdaMax(params::ParamType{Tdim},
                      qL::AbstractVector{Tsol}) where {Tsol, Tdim}

  gamma = params.gamma
  Un = zero(Tsol)
  rhoLinv = 1/qL[1]

  pL = calcPressure(params, qL)
  aL = sqrt(gamma*pL*rhoLinv)  # speed of sound

  for i=1:Tdim
    u_i = qL[i+1]*rhoLinv
    Un += u_i*u_i
  end

  lambda_max = sqrt(Un) + aL

  return lambda_max
end

"""
  Compute approximate wave speeds for a Riemann problem using the 2 rarefaction
  model (see Toro's Riemann Solvers book).
  These values are upper bounds for the true wave speeds for 1 <= gamma <= 5/3.

  Note that the wave speeds are *not* scaled by the length of the normal
  vector.

  **Inputs**

   * params
   * qL
   * qR
   * nrm: normal vector

  **Outputs**

   * sL: slowest wave speed
   * sR: fastest wave speed
"""
function calc2RWaveSpeeds(params::ParamType{Tdim}, qL::AbstractVector{Tsol},
                          qR::AbstractVector, nrm::AbstractVector{Tmsh}
                         ) where {Tdim, Tsol, Tmsh}

  Tres = promote_type(Tsol, Tmsh)

  pL = calcPressure(params, qL); pR = calcPressure(params, qR)
  aL = sqrt(params.gamma*pL/qL[1]); aR = sqrt(params.gamma*pR/qR[1])
  z = params.gamma_1/(2*params.gamma)

  # compute velocity in face normal direction
  u_nrmL = zero(Tres); u_nrmR = zero(Tres)
  fac = calcLength(params, nrm)
  for i=1:Tdim
    u_nrmL += qL[i+1]*nrm[i]
    u_nrmR += qR[i+1]*nrm[i]
  end
  u_nrmL /= fac*qL[1]; u_nrmR /= fac*qR[1]

  num = aL + aR - 0.5*params.gamma_1*(u_nrmR - u_nrmL)
  den = aL/(pL^z) + aR/(pR^z)

  p_tr = (num/den)^(1/z)

  # compute q values
  if p_tr <= pL
    qfL = Tres(1.0)
  else
    qfL = sqrt(1 + (params.gamma + 1)*(p_tr/pL - 1)/(2*params.gamma))
  end

  if p_tr <= pR
    qfR = Tres(1.0)
  else
    qfR = sqrt(1 + (params.gamma + 1)*(p_tr/pR - 1)/(2*params.gamma))
  end

#  qfL = 1.3
#  qfR = 1.02

  sL = u_nrmL - aL*qfL
  sR = u_nrmR + aR*qfR

  return sL, sR
end

function calcMomentContribution!(sbpface::AbstractFace{Tsbp}, xsbp::AbstractArray{Tmsh,3},
    dforce::AbstractArray{Tres,3}, xyz_about::AbstractArray{Tmsh,1}) where {Tsbp,Tmsh,Tres
                                             }
  @assert( sbpface.numnodes == size(xsbp,2) == size(dforce,2) )
  @assert( size(xsbp,3) == size(dforce,3) )
  @assert( size(xsbp,1) == size(dforce,1) == size(xyz_about,1) )
  rvec = zeros(Tmsh,3)
  dM = zeros(Tres,3)
  moment = zeros(Tres, length(xyz_about))
  for f = 1:size(dforce,3)
    for i = 1:sbpface.numnodes
      for di = 1:3
        rvec[di] = xsbp[di,i,f] - xyz_about[di]
      end
      crossProd(rvec, sview(dforce,:,i,f), dM)
      for di = 1:3
        moment[di] += dM[di]*sbpface.wface[i]
      end
    end
  end

  return moment
end

function calcMomentContribution!(mesh::AbstractMesh{Tmsh},
               eqn::AbstractSolutionData{Tsol, Tres},
               bndry_nums::Array{Int, 1},
               xyz_about::AbstractArray{Tmsh, 1}) where {Tmsh, Tsol, Tres}

  moment = zeros(Tres, mesh.dim)
  for i=1:length(bndry_nums)
    start_idx = mesh.bndry_offsets[ bndry_nums[i] ]
    end_idx = mesh.bndry_offsets[ bndry_nums[i] + 1 ] - 1
    face_range = start_idx:end_idx
    bndry_faces = sview(mesh.bndryfaces, face_range)
    coords = ro_sview(mesh.coords_bndry, :, :, face_range)
    nrm = ro_sview(mesh.nrm_bndry, :, :, face_range)

    # compute dforce
    dforce = computeDForce(mesh, eqn, bndry_faces, nrm)

    # compute moment
    moment += calcMomentContribution!(mesh.sbpface, coords, dforce, xyz_about)

  end

  return moment
end

function calcMomentContribution_revm!(mesh::AbstractMesh, eqn::AbstractSolutionData, bndry_nums::Array{Int, 1}, xyz_about::AbstractArray{Tmsh, 1}, moment_bar::AbstractArray{Tres, 1}) where {Tmsh, Tres}

  nfaces = length(mesh.bndryfaces)
  moment = zeros(Tres, mesh.dim)
  for i=1:length(bndry_nums)
    start_idx = mesh.bndry_offsets[ bndry_nums[i] ]
    end_idx = mesh.bndry_offsets[ bndry_nums[i] + 1 ] - 1
    face_range = start_idx:end_idx
    bndry_faces = ro_sview(mesh.bndryfaces, face_range)
    coords = ro_sview(mesh.coords_bndry, :, :, face_range)
    coords_bar = zeros(coords)

    nrm = sview(mesh.nrm_bndry, :, :, face_range)
    nrm_bar = sview(mesh.nrm_bndry_bar, :, :, face_range)

    # compute dforce
    dforce = computeDForce(mesh, eqn, bndry_faces, nrm)
    dforce_bar = zeros(dforce)

    #--------------------------------------------------------------------------
    # start reverse sweep
    calcMomentContribution_rev!(mesh.sbpface, coords, coords_bar, dforce, dforce_bar, xyz_about, moment_bar)

    computeDForce_revm!(mesh, eqn, bndry_faces, nrm_bar, dforce_bar)
  end

  return nothing
end

function computeDForce(mesh::AbstractMesh,
                       eqn::AbstractSolutionData{Tsol, Tres},
                       bndryfaces::AbstractArray{Boundary, 1},
                       nrm::Abstract3DArray{Tmsh}) where {Tmsh, Tsol, Tres}

  nfaces = length(bndryfaces)
  dforce = zeros(Tres, mesh.dim, mesh.numNodesPerFace, nfaces)
  for i=1:nfaces
    for j=1:mesh.numNodesPerFace
      q_j = ro_sview(eqn.q_bndry, :, j, i)  # q_bndry must already have been populated
      p = calcPressure(eqn.params, q_j)
      for k=1:mesh.dim
        dforce[k, j, i] = p*nrm[k, j, i]
      end
    end
  end

  return dforce
end

function computeDForce_revm!(mesh::AbstractMesh{Tmsh},
           eqn::AbstractSolutionData{Tsol, Tres},
           bndryfaces::AbstractArray{Boundary, 1},
           nrm_bar::Abstract3DArray,
           dforce_bar::Abstract3DArray) where {Tmsh, Tsol, Tres}

   nfaces = length(bndryfaces)
   for i=1:nfaces
    for j=1:mesh.numNodesPerFace
      q_j = ro_sview(eqn.q_bndry, :, j, i)
      p = calcPressure(eqn.params, q_j)
      for k=1:mesh.dim
        nrm_bar[k, j, i] = dforce_bar[k, j, i]*p
      end
    end
  end

  return nothing
end

@doc """
### calcMomentContribution_rev!

This is the reverse differentiated version of calcMomentContribution!.  See docs
of calcMomentContribution! for further details of the primal method.  This
function is differentiated with respect to the primal version's `xsbp` and
`dforce` variables.

**Inputs**

* `sbpface`: an SBP face operator type
* `xsbp`: SBP-face nodes in physical space; [coord, sbp node, face]
* `dforce`: scaled force at the sbpface nodes; [coord, sbp node, face]
* `xyz_about`: point about which the moment is taken
* `moment_bar`: left multiplies d(moment)/d(xsbp) and d(moment)/d(dforce)

**InOuts**

* `xsbp_bar`: result of vector Jacobian product; [coord, sbp node, face]
* `dforce_bar`: result of vector Jacobian product; [coord, sbp node, face]

"""->
function calcMomentContribution_rev!(sbpface::AbstractFace{Tsbp}, xsbp::AbstractArray{Tmsh,3},
    xsbp_bar::AbstractArray{Tmsh,3}, dforce::AbstractArray{Tsol,3},
    dforce_bar::AbstractArray{Tsol,3}, xyz_about::AbstractArray{Tmsh,1},
    moment_bar::AbstractArray{Tres,1}) where {Tsbp,Tmsh,Tsol,Tres
           }
  @assert( sbpface.numnodes == size(xsbp,2) == size(xsbp_bar,2)
           == size(dforce,2) == size(dforce_bar,2) )
  @assert( size(xsbp,3) == size(xsbp,3) == size(dforce,3) == size(dforce_bar,3) )
  @assert( size(xsbp,1) == size(xsbp_bar,1) == size(dforce,1)
           == size(dforce_bar,1) == size(xyz_about,1) )
  rvec = zeros(Tmsh,3)
  rvec_bar = zeros(Tres,3)
  dM_bar = zeros(Tres,3)
  for f = 1:size(dforce,3)
    for i = 1:sbpface.numnodes
      for di = 1:3
        rvec[di] = xsbp[di,i,f] - xyz_about[di]
        rvec_bar[di] = 0
        dM_bar[di] = 0  # ???
      end
      # start reverse sweep
      for di = 1:3
        # moment[di] += dM[di]*sbpface.wface[i]
        dM_bar[di] += moment_bar[di]*sbpface.wface[i]
      end
      # crossProd(rvec, sview(dforce,:,i,f), dM)
      crossProd_rev(rvec, rvec_bar, sview(dforce,:,i,f),
                       sview(dforce_bar,:,i,f), dM_bar)
      for di = 1:3
        # rvec[di] = xsbp[di,i,f] - xyz_about[di]
        xsbp_bar[di,i,f] += rvec_bar[di]
      end
    end
  end
end

#=
function calcMomentContribution!(sbpface::AbstractFace{Tsbp},
    xsbp::AbstractArray{Tmsh,3},
    dforce::AbstractArray{Tres,3}, xyz_about::AbstractArray{Tmsh,1}) where {Tsbp,Tmsh,Tres}
  @assert( sbpface.numnodes == size(xsbp,2) == size(dforce,2) )
  @assert( size(xsbp,3) == size(dforce,3) )
  @assert( size(xsbp,1) == size(dforce,1) == size(xyz_about,1) )
  rvec = zeros(Tmsh,3)
  dM = zeros(Tres,3)
  moment = zeros(Tres, length(xyz_about))
  for f = 1:size(dforce,3)
    for i = 1:sbpface.numnodes
      for di = 1:3
        rvec[di] = xsbp[di,i,f] - xyz_about[di]
      end
      crossProd(rvec, sview(dforce,:,i,f), dM)
      for di = 1:3
        moment[di] += dM[di]*sbpface.wface[i]
      end
    end
  end

  return moment
end
=#


text:
@testset "Closeness" begin
    g5 = DiGraph(4)
    add_edge!(g5,1,2); add_edge!(g5,2,3); add_edge!(g5,1,3); add_edge!(g5,3,4)
    for g in testdigraphs(g5)
      y = @inferred(closeness_centrality(g; normalize=false))
      z = @inferred(closeness_centrality(g))
      @test y == [0.75, 0.6666666666666666, 1.0, 0.0]
      @test z == [0.75, 0.4444444444444444, 0.3333333333333333, 0.0]
    end

    g5 = Graph(5)
    add_edge!(g5,1,2)
    for g in testgraphs(g5)
      z = @inferred(closeness_centrality(g))
      @test z[1] == z[2] == 0.25
      @test z[3] == z[4] == z[5] == 0.0
    end
end

ext: jl

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
# functions to preprocess the data and classify some information of model that algorithm needed
include("modelPreprocessing.jl")
# functions of algorithm
include("algorithm.jl")

# solveDEA: solving dea problems by a smart algorithm
function solveDEA(m::JuMP.Model; lpUB=300)
  # parameters
  # input
  # m: model variable of JuMP which is entered by JuMP's scripts
  # lpUB: the limit size of LP

  # data preprocessing
  # getting the data of full-size model which is entered by users
  fullsizeObjSense, fullsizeObjCoeff, fullsizeConstLB, fullsizeConstUB, fullsizeConstMatrix, fullsizeVarLB, fullsizeVarUB, fullsizeNumVar, fullsizeNumConst = getFullsizeModelData(m)

  # checking the formulation of model is envelopment form or mutiplier form
  # envelopment form and mutiplier form are dual relationships
  modelForm = identifyModelForm(fullsizeObjCoeff)

  # if model is mutiplier form, the model has to transform to envelopment form
  if modelForm == :Multiplier
    # tranforming model to envelopment form
    fullsizeObjSense, fullsizeObjCoeff, fullsizeConstLB, fullsizeConstUB, fullsizeConstMatrix, fullsizeVarLB, fullsizeVarUB, fullsizeNumVar, fullsizeNumConst = modelTransform(fullsizeObjSense, fullsizeObjCoeff, fullsizeConstLB, fullsizeConstUB, fullsizeConstMatrix, fullsizeVarLB, fullsizeVarUB, fullsizeNumVar, fullsizeNumConst)
  end

  # finding the index of theta w.r.t the variable vector
  thetaIdx = findThetaIdx(fullsizeObjCoeff)
  # getting the target DMU by rhs and the coefficients of theta in constraints
  targetDMU = getTargetDMU(fullsizeConstLB, fullsizeConstUB, fullsizeConstMatrix, thetaIdx)
  # finding the index of target DMU w.r.t the variable vector
  targetDMUIdx = findTargetDMUIdx(fullsizeConstMatrix, targetDMU)
  # end data preprocessing

  # randomly selecting DMUs to build a small model with predefined lp size of lpUB
  set, notset, smallConstMatrix = sampling(lpUB, fullsizeConstMatrix, thetaIdx, targetDMUIdx)

  #generating the data of small model by set
  smallVarLB = generateSmallV(fullsizeVarLB, lpUB, set, targetDMUIdx, thetaIdx)
  smallVarUB = generateSmallV(fullsizeVarUB, lpUB, set, targetDMUIdx, thetaIdx)
  smallObjCoeff = generateSmallV(fullsizeObjCoeff, lpUB, set, targetDMUIdx, thetaIdx)
  smallConstLB = fullsizeConstLB
  smallConstUB = fullsizeConstUB
  smallSense = fullsizeObjSense

  # definition and initialization
  optStatus = false   # the optimality status of model
  iterations = 0      # the iterations that algorithm uses to solve model
  overNum = 0         # the most infeasible DMU, which is the indices w.r.t. notset

  lm = MathProgBase.LinearQuadraticModel(m.solver)  # model variable of MathProgBase, we use this to build small model
  stat = :Default     # model status
  m.objVal = NaN      # objective value of full-size model
  m.colVal = fill(NaN, fullsizeNumVar)    # values of all variables of full-size model
  m.linconstrDuals = Array(Float64, 0)    # dual values of full-size model

  # loop until model being solved to optimal or algorithm processing over 1000 iterations
  while optStatus == false
    iterations = iterations + 1
    if iterations > 1000
      warn("Iterations > 1000")
      stat = :Unsolved
      break
    end
    optStatus = true

    lm = MathProgBase.LinearQuadraticModel(m.solver)
    # loading model information
    MathProgBase.loadproblem!(lm, smallConstMatrix, smallVarLB, smallVarUB, smallObjCoeff, smallConstLB, smallConstUB, smallSense)
    # solving model
    MathProgBase.optimize!(lm)

    stat = MathProgBase.status(lm)
    # getting dual values of small model
    constDuals = MathProgBase.getconstrduals(lm)

    # checking optimality and getting the most infeasible DMU which will be added to set
    optStatus, overNum = checkOpt(constDuals, notset, fullsizeConstMatrix, smallSense)

    # redefining set if the model is not optimal
    if optStatus == false
      # adding the most infeasible DMU to set and deleting a DMU from set by the randomly selecting order
      # updating model
      reSampling(overNum, iterations, set, notset, smallConstMatrix, smallVarLB, smallVarUB, smallObjCoeff,
                 fullsizeConstMatrix, fullsizeVarLB, fullsizeVarUB, fullsizeObjCoeff)
    end
  end

  # Extracting solution
  if stat == :Optimal
    # recording objective value
    m.objVal = MathProgBase.getobjval(lm) + m.obj.aff.constant

    if modelForm == :Multiplier
      m.linconstrDuals = fill(0.0, fullsizeNumVar)
      solution = MathProgBase.getsolution(lm)

      # recording dual value
      m.linconstrDuals[set] = solution[1:length(set)]
      m.linconstrDuals[targetDMUIdx] = solution[length(set)+1]
      m.linconstrDuals[thetaIdx] = solution[length(set)+2]

      # recording value of variables
      m.colVal = try
          MathProgBase.getconstrduals(lm)[1:fullsizeNumConst]
      catch
          fill(NaN, fullsizeNumConst)
      end
    else
      m.colVal = fill(0.0, fullsizeNumVar)
      solution = MathProgBase.getsolution(lm)

      # recording value of variables
      m.colVal[set] = solution[1:length(set)]
      m.colVal[targetDMUIdx] = solution[length(set)+1]
      m.colVal[thetaIdx] = solution[length(set)+2]

      # recording dual value
      m.linconstrDuals = try
          MathProgBase.getconstrduals(lm)[1:fullsizeNumConst]
      catch
          fill(NaN, fullsizeNumConst)
      end
    end
  end

  return stat, iterations
end


text:
include("timeseries.jl")
include("datasets.jl")
abstract type DataGroup end

mutable struct SeriesDataGroup <:DataGroup #signal равномерно дискретный ряд
    filepath::String
    groupname::String

    TG::TimeGrid
    A::PhysicalInfo
    data
    result #итоговый результат всех правок над данными
    UndoRedo::UndoRedoTool
end

mutable struct EventDataGroup <:DataGroup #события
    filepath::String
    groupname::String

    TG::TimeGrid
    indexdata::String
    data
    result #итоговый результат всех правок над данными
    UndoRedo::UndoRedoTool
end

mutable struct ChannelDataGroup <:DataGroup # канал (все данные, рассчитанные по приборному каналу)
    filepath::String
    groupname::String

    TG::TimeGrid
    A::PhysicalInfo
    data #::StructArray
    result #итоговый результат всех правок над данными
    UndoRedo::UndoRedoTool
end

mutable struct SegmentDataGroup <:DataGroup #сегмент
    # filepath::String
    groupname::String

    TG::TimeGrid
    ibegdata:: String
    ienddata:: String
    typename:: String
    mask:: Dict

    data
    result #итоговый результат всех правок над данными
    UndoRedo:: UndoRedoTool

end

mutable struct UnknownDataGroup <:DataGroup #сегмент
    filepath::String
    groupname::String

    attr::TimeGrid

    data
    result #итоговый результат всех правок над данными
    UndoRedo::UndoRedoTool

end

function dg_new(baseIP::IPv4,port::Union{String,Int64},groupName::String)

    r = HTTP.request("GET", "http://$baseIP:$port/api/getDataTree")
    tree = JSON.parse(String(r.body))[1]
    node =  findnode(tree, groupName) #запрос группы

    attr = node["attrs"]

    #пока берем так. но надо бы из группы данных
    r = HTTP.request("GET", "http://$baseIP:$port/api/getData?dataName=Freq&index=0&from=0&count=1")
    Freq = reinterpret(Int32, base64decode(r.body)) |> collect
    r = HTTP.request("GET", "http://$baseIP:$port/api/getData?dataName=StartTime&index=0&from=0&count=1")
    StartTime = reinterpret(Int32, base64decode(r.body)) |> collect

    attr["Freq"] = Freq[1]
    attr["StartTime"] = StartTime[1]
    #определяем временную сетку. хз зачем, правда
    TG,PhInfo = parseAttr(attr)

    #создаем пустую штуку для хранения истории
    URT = UndoRedoTool()
    #если это не группа(а вдруг датасет!)
    if !haskey(attr,"grouptype")
        return []
    else
        if haskey(attr, "grouptype") && attr["grouptype"]=="segment"
            ibegdata = attr["ibegdata"]
            ienddata = attr["ienddata"]
            if haskey(attr,"typedata")
                featureName = attr["typedata"]
                loadDsNames = [ibegdata, ienddata, featureName]
            else
                featureName = "none"
                loadDsNames = [ibegdata, ienddata]
            end

            mask = getMask(parseType(attr))
            if isempty(mask)
                #если не было маски, то называем тип по имени датагруппы
                mask[groupName] = AllMask(UInt32(1),UInt32(1))
            end

            DS = Dict{String,Any}()
            for ds in loadDsNames
                node =  findnode(tree, ds) #запрос группы
                attrds = node["attrs"]
                attrds["datatype"] = eval(Symbol(node["type"]))
                attrds["Freq"] = Freq[1]
                attrds["StartTime"] = StartTime[1]
                DS[ds] = ds_new(baseIP,port,groupName,ds, attrds)
            end
            # DataSets, data = loadDataSets(fid,datapath)
            # dataStr = StructArray((name = DataSets, data = data))
            # if any(isa.(data,FeatureDataSet))
            #     featureName = DataSets[isa.(data,FeatureDataSet)][1]
            # else
            #     featureName = "none"
            # end
             DG_obj = SegmentDataGroup(groupName,TG,ibegdata,ienddata,featureName,mask,DS,[],URT)
        # elseif attr["grouptype"]=="series"
        #
        #     DataSets, data = loadDataSets(fid,datapath)
        #     dataStr = StructArray((name = DataSets, data = data))
        #
        #     DG_obj = SeriesDataGroup(filepath, datapath, fid, TG, PhInfo, dataStr,[], URT)
        # elseif attr["grouptype"]=="event"
        #     index = attr["indexdata"]
        #
        #     DataSets, data = loadDataSets(fid,datapath)
        #     dataStr = StructArray((name = DataSets, data = data))
        #
        #     DG_obj = EventDataGroup(filepath, datapath, fid, TG, index, dataStr,[], URT)
        #
        # elseif attr["grouptype"]=="channel"
        #
        #     DataSets, data = loadDataSets(fid,datapath)
        #     dataStr = StructArray((name = DataSets, data = data))
        #     DG_obj = ChannelDataGroup(filepath, datapath, fid, TG, PhInfo, dataStr,[], URT)

        else #непонятные данные
            # DataSets, data = loadDataSets(fid,datapath)
            # dataStr = StructArray((name = DataSets, data = data))
            # DG_obj = UnknownDataGroup(filepath, datapath, fid, attr, dataStr, [], URT)
            DG_obj = []
        end
    end
    DG_obj.UndoRedo.sourse = DG_obj
    return DG_obj
end

#смотрим, исходные ширины это сегменты с началом-концом или шириной
function getSegBegsType(DG::SegmentDataGroup)

    ibeg =  DG.data[DG.ibegdata].data[1]("0","","")
    iendDS = DG.data[DG.ienddata]
    # @info size(ibeg)
    typename = DG.typename

    isW = isa(iendDS, IntervalDataSet) #ширина ли это
    if isW
        iend = ibeg + Int32.(iendDS.data[1]("0","",""))
    else
        iend = iendDS.data[1]("0","","")
    end
    if occursin(typename,"none") #none там, где нет данных типа
        type = ones(Int64,size(ibeg))
    else
        type =DG.data[typename].data[1]("0","","")
    end

    ibeg, iend, type, isW
end
#добавление нового сегмента
#command тут не обрабатывается, только пишется в историю правок датагруппы
function addSeg!(DG::SegmentDataGroup,newSeg::StructArray, mode::String)
    if isempty(DG.result) || DG.UndoRedo.state==0 #&& (!haskey(URT.result["ibeg"]) || isempty(URT.result["ibeg"]))
        #если ничего не делалось над объектом, то читаем данные из источника
        ibeg, iend, type, isW = getSegBegsType(DG)
        @info sizeof(ibeg), sizeof(iend), sizeof(type)
        oldSeg = StructArray(ibeg = ibeg, iend = iend, type = type)
    else
        #если были правки раньше
        iendDS =  DG.data[DG.ienddata]
        isW = isa(iendDS, IntervalDataSet) #ширина ли это

        resultOld = DG.result

        list = getfield(resultOld, :fieldarrays) #тк в резалте другие имена, пересобираем в стандартные
        if isW
            iend = list[1] + list[2]
        else
            iend = list[2]
        end
        oldSeg = StructArray((list[1], iend, list[3]), names = (Symbol("ibeg"),Symbol("iend"),Symbol("type")))
    end
    result = add_seg(oldSeg,newSeg,mode)
    # @info result
    if isW
        iend = result.iend - result.ibeg
    else
        iend = result.iend
    end
    println("пишем резаль")
    DG.result = StructArray((convert.(DG.data[DG.ibegdata].datatype,result.ibeg), convert.(DG.data[DG.ienddata].datatype,iend), result.type), names = (Symbol(DG.ibegdata),Symbol(DG.ienddata),Symbol(DG.typename)))
    # @info DG.result
end

#изменение типа сегментов, попавших в диапазон от-до
function changeType!(DG::SegmentDataGroup,from::Int64,to::Int64, mode::String,type_names::Union{String,Vector{String}},mask::Dict)
    if DG.UndoRedo.state==0 #&& (!haskey(URT.result["ibeg"]) || isempty(URT.result["ibeg"]))
        #если ничего не делалось над объектом, то читаем данные из источника
        #если ничего не делалось над объектом, то читаем данные из источника
        ibeg, iend, type, isW = getSegBegsType(DG)
    else
        #если были правки раньше
        iendDS =  DG.data[DG.ienddata]
        isW = isa(iendDS, IntervalDataSet) #ширина ли это

        resultOld = DG.result
        list = getfield(resultOld, :fieldarrays) #тк в резалте другие имена, пересобираем в стандартные
        if isW
            iend = list[1] + list[2]
        else
            iend = list[2]
        end
        ibeg = list[1]
        type = list[3]
    end
    #определяем попавшие в диапазон сегменты
    segRange = searchinrange(ibeg,iend,from,to)
    if mode == "delete"
        #прям удаляем сегменты, попавшие под выделение
        deleteat!(ibeg,segRange)
        deleteat!(iend,segRange)
        deleteat!(type,segRange)
    elseif mode=="setbit"
        #присваеваем новый тип
        new_types = setbit(mask,type[segRange],type_names)
        type[segRange] .= new_types
    elseif mode=="rewrite"
        new_types = detIntType(DG.mask,type_names)
        type[segRange] .= new_types
    end
    if isW
        iend = iend - ibeg
    end
    #записываем новые сегменты
    DG.result = StructArray((convert.(DG.data[DG.ibegdata].datatype,ibeg), convert.(DG.data[DG.ibegdata].datatype,iend), type), names = (Symbol(DG.ibegdata),Symbol(DG.ienddata),Symbol(DG.typename)))
end

function parseCommand(baseIP::IPv4,port::Union{String,Int64},DG::DataGroup,command::Vector{String},flUR = 0)
    #можно подавать команды в векторе - будут применены последовательно
    for com in command
        parseCommand(baseIP, port, DG,com,flUR)
    end
end
function parseCommand(baseIP::IPv4,port::Union{String,Int64},AllObj::Dict,command::Vector{String},flUR = 0)
    #можно подавать команды в векторе - будут применены последовательно
    for com in command
        parseCommand(baseIP, port, AllObj,com,flUR)
    end
end

function parseCommand(baseIP::IPv4,port::Union{String,Int64},DG::DataGroup,command::String,flUR = 0)
    #проверяем, над этой ли датагруппой команда
    manualEvent = JSON.parse(command)
    #определяем имя данных
    gp_name = manualEvent["chName"]
    tempDict = Dict{String,Any}()
    tempDict["history"] = []
    tempDict["state"] = 0
    tempDict["dataStorage"] = Dict{String,Any}()
    tempDict["dataStorage"][gp_name] = DG

    parseCommand(baseIP, port, tempDict,command,flUR)
end

#обработка команды-правки пользователя
#AllObj - набор датагрупп в памяти
#command - прочитанный в String JSON-file
#flUR - флаг Undo-Redo , если true, то история не перетирается
function parseCommand(baseIP::IPv4,port::Union{String,Int64}, AllObj::Dict,command::String,flUR = 0)
    #разбор JSON-a
    #судя по докам, он умеет парсить только в словарь
    #на сервере уже разбирали файл-команду
    manualEvent = JSON.parse(command)
    #пока считаем, что имя группы лежит в chName
    gp_name = manualEvent["chName"]
    # @info gp_name
    comandID = manualEvent["command"]["id"] #определяем, какую команду делаем с сегментом
    #в sessionID идентификатор события
    if haskey(manualEvent,"sessionID")
        sessionID = manualEvent["sessionID"]
    else
        sessionID = "12345"
    end
    #ТУТ НАДО НАЙТИ ОБЪЕКТ, КОТОРЫЙ СООТВЕТСТВУЕТ ДАННЫМ, КОТОРЫЕ ПРАВИМ
    if ~isempty(AllObj["dataStorage"]) && haskey(AllObj["dataStorage"],gp_name)
        DG = AllObj["dataStorage"][gp_name]
    else
        DG = dg_new(baseIP,port,gp_name)
        # @info AllObj["dataStorage"]
        AllObj["dataStorage"][gp_name] = DG
        #тут можно использовать targetData
    end
    # println(DG)

    # filePath = collect(splitpath(DG.filepath))

    #Это бы надо куда-то снаружи, но ладно. Мы ж независимы..
    filePath =  sessionID*"_history.json"
    historyPath = sessionID*"_history.json"
    # @info historyPath
    if comandID == "ADD_SEGMENT"
        #границы добавляемого фрагмента
        ibeg = convert.(Int,manualEvent["command"]["args"]["ibeg"])
        iend = convert.(Int,manualEvent["command"]["args"]["iend"])
        type = convert.(String,manualEvent["command"]["args"]["type"])
        if type == "+"
            type = string(keys(DG.mask))[3:end-2]
        end
        if haskey(manualEvent["command"]["args"],"mode")
            mode = manualEvent["command"]["args"]["mode"]
        else
            mode = "rewrite"
        end
        typeNum = detIntType(DG.mask,type) #преобразуем тип в число
        type = Vector{Int}(undef,length(ibeg)) #если несколько
        type[1:end] .= typeNum

        if !isa(ibeg,Vector) #если было одно значение только, то его в вектор надо
            ibeg = [ibeg]
            iend = [iend]
        end
        println("Добавляем сегмент: ")
        @info ibeg, iend, type
        newSeg = StructArray(ibeg = ibeg, iend = iend, type = type)
        #addSeg!(DG,newSeg,"rewrite",commandToString(manualEvent["command"]))
        addSeg!(DG,newSeg,mode)

    elseif comandID == "CHANGE_TYPE" || comandID == "DELETE_SEGMENT"
        from = manualEvent["command"]["args"]["ibeg"][1]
        to = manualEvent["command"]["args"]["iend"][1]
        newType = manualEvent["command"]["args"]["type"]
        println("Удаляем/меняем тип сегмент: ")

        @info from, to

        if comandID == "CHANGE_TYPE"
            mode = manualEvent["command"]["args"]["mode"]
        else
            mode = "delete"
            newType = Vector{String}()
        end
        changeType!(DG,from,to,mode,newType,DG.mask) #type - вектор String c именами
    else
        println("Unknown command")
    end
    #обновляем состояние,наращиваем историю
    #если стейт и количество объектов в истории совпадают - то все норм
    if flUR==0
        #это не  undo-redo - пришла новая правка
        if DG.UndoRedo.state==size(DG.UndoRedo.history,1)
            DG.UndoRedo.state+=1
            push!(DG.UndoRedo.history,command)

            push!(AllObj["history"],command)
            AllObj["state"] += 1
        elseif DG.UndoRedo.state<=size(DG.UndoRedo.history,1)
            #если была новая правка после undo, То все после в истории надо перетереть
            deleteat!(DG.UndoRedo.history,DG.UndoRedo.state+1:size(DG.UndoRedo.history,1))
            DG.UndoRedo.state+=1
            push!(DG.UndoRedo.history,command)

            deleteat!(AllObj["history"],AllObj["state"]+1:size(AllObj["history"],1))
            push!(AllObj["history"],command)
            AllObj["state"] += 1
        end
        # пишем
        writeManualMark(command, historyPath)

    elseif flUR == 1
        #1 - значит redo
        # пишем команду
        writeManualMark(command, historyPath)

    elseif flUR == -1
        #-1 - значит undo и надо удалить последнюю строку из файла
        #ничего не делаем
        deleteLastCommand(historyPath)
    end
    #если были заданы границы, то выдаем сегменты в них
    if haskey(manualEvent,"fromto")
        segInRange = getSegInRange(DG, manualEvent["fromto"][1],manualEvent["fromto"][2])
        return segInRange
    else
        return AllObj
    end

end

function findSegDS(DG::SegmentDataGroup,from::Union{Int32,Int64},to::Union{Int32,Int64},features="")

    if !isempty(DG.result) && DG.UndoRedo.state>0
        iendDS =  DG.data[DG.ienddata]
        isW = isa(iendDS, IntervalDataSet) #ширина ли это

        resultOld = DG.result
        list = getfield(resultOld, :fieldarrays) #тк в резалте другие имена, пересобираем в стандартные
        if isW
            iend = list[1] + list[2]
        else
            iend = list[2]
        end
        ibeg = list[1]
        type = list[3]

        segs = DG.result
    else
        ibeg, iend, type, isW = getSegBegsType(DG)
        if isW
            iendW = iend - ibeg
        else
            iendW = iend
        end
        segs = StructArray(ibeg = ibeg,iend = iendW,type = type)
    end

    return segs,ibeg,iend
end
#запрос сегментов в диапазоне - нужно для отправки клиенту для отрисовки
#в features добавить фильтр по типу
function getSegInRange(DG::SegmentDataGroup,from::Int64,to::Int64,features="")

    segs,ibeg,iend = findSegDS(DG,from,to)
    range = searchinrange(ibeg,iend,from,to)
    SIR = segs[range]
    """
    добавить фильтр по типам из features
    """
    #тут надо наверное делать json или еще что-то
    return SIR
end

function getStructData(DG::Union{SegmentDataGroup,EventDataGroup},from,to,dsTake="all")
    if isa(DG,SegmentDataGroup)
        segs,ibeg,iend = findSegDS(DG,from,to)
        range = searchinrange(ibeg,iend,from,to)
        data = getData(DG,range.start,range.stop,dsTake)

        return data
    end
end
#запрос данных из датагрупп
function getData(DG::Union{SegmentDataGroup,EventDataGroup},from,to,dsTake="all")
    ind = Int32(from):Int32(to)
    allInRes = true
    data = Vector{Any}()
    dsTake = split(dsTake,',')
    if dsTake != "all"
        if ~isa(dsTake,Vector)
            dsTake = [dsTake]
        end
        for dsT in dsTake
            #если хотим что-то из нередактируемых данных
            allInRes = allInRes && any(Symbol(dsT).==propertynames(DG.result)) && dsT!="original"
        end
    end
    if !isempty(DG.result) && allInRes
        allDS = DG.result

        list = getfield(allDS, :fieldarrays)
        if dsTake=="all"
            for i = 1:length(propertynames(allDS))
                addData = list[i][ind]
                push!(data,addData)
            end
        else
            i=1
            for key in dsTake
                fi = findall(x-> x== Symbol(key), propertynames(allDS))
                if !isempty(fi) #запрашеваемое поле найдено
                    addData = convert.(DG.data[key].datatype,list[fi[1]][ind])
                    # addData = list[fi[1]][ind]
                    push!(data,addData)
                end
                i+=1
            end
            if length(dsTake) == 1
                data = data[1]
            end
        end
        return data
    else
        ##ДОБАВИТЬ ЧТЕНИЕ ДАННЫХ ИЗ ИСХОДНИКА
        if dsTake=="all"
            for key in keys(DG.data)
                addData = DG.data[key].data[1](string(from),srting(to),"")
                push!(data,addData)
            end
        else
            i=1
            for key in dsTake
                if haskey(DG.data, key)
                    addData = DG.data[key].data[1](string(from),string(to+1),"")
                    push!(data,addData)
                end
                i+=1
            end
        end

        return data

        # allDS = DG.data.data
        # DSs = [[DG.data.data[i].data] for i = 1:length(DG.data)]
        # allDS = StructArray((DSs...,), names = (Symbol.(DG.data.name)...,))
    end

end

#отмена для датагруппы
function undo!(DG::DataGroup,baseIP::IPv4,port::Union{String,Int64},)
    #уперлись в левый край - ничего не делаем
    if DG.UndoRedo.state == 0
        return DG
    end
    #понижаем статус
    DG.UndoRedo.state -= 1
    #очищаем результаты
    DG.result = []
    if DG.UndoRedo.state != 1 && DG.UndoRedo.state != 0
        #повторяем все команды до стейта с флагом "undo-redo"
        for command in DG.UndoRedo.history[1:DG.UndoRedo.state-1]
            parseCommand(baseIP, port, DG,command,-2)
        end
        #только последнюю отдельно, чтобы удалить не забыть строчку в файле правок
        parseCommand(baseIP, port, DG,DG.UndoRedo.history[DG.UndoRedo.state],-1)
    elseif DG.UndoRedo.state == 1
        parseCommand(baseIP, port, DG, DG.UndoRedo.history[DG.UndoRedo.state],-1)
    end
    return DG
end

#возврат для датагруппы
function redo!(DG::DataGroup,baseIP::IPv4,port::Union{String,Int64})
    #уперлись в правый край - ничего не делаем
    if DG.UndoRedo.state == size(DG.UndoRedo.history,1)
        return DG
    end
    #повышаем статус
    DG.UndoRedo.state += 1
    #повторяем команду этого стейта
    command = DG.UndoRedo.history[DG.UndoRedo.state]
    parseCommand(baseIP, port, DG,command,1)
    return DG
end


text:
brusselator_f(x, y, t) = ifelse((((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) &&
                                (t >= 1.1), 5., 0.)
function limit(a, N)
  if a == N+1
    return 1
  elseif a == 0
    return N
  else
    return a
  end
end
function brusselator_2d_loop(du, u, p, t)
  @inbounds begin
    A, B, α, xyd, dx, N = p
    α = α/dx^2
    for I in CartesianIndices((N, N))
      x = xyd[I[1]]
      y = xyd[I[2]]
      i = I[1]
      j = I[2]
      ip1 = limit(i+1, N)
      im1 = limit(i-1, N)
      jp1 = limit(j+1, N)
      jm1 = limit(j-1, N)
      du[i,j,1] = α*(u[im1,j,1] + u[ip1,j,1] + u[i,jp1,1] + u[i,jm1,1] - 4u[i,j,1]) +
      B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y, t)
    end
    for I in CartesianIndices((N, N))
      i = I[1]
      j = I[2]
      ip1 = limit(i+1, N)
      im1 = limit(i-1, N)
      jp1 = limit(j+1, N)
      jm1 = limit(j-1, N)
      du[i,j,2] = α*(u[im1,j,2] + u[ip1,j,2] + u[i,jp1,2] + u[i,jm1,2] - 4u[i,j,2]) +
      A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]
    end
  end
end
function init_brusselator_2d(xyd)
  N = length(xyd)
  u = zeros(N, N, 2)
  for I in CartesianIndices((N, N))
    x = xyd[I[1]]
    y = xyd[I[2]]
    u[I,1] = 22*(y*(1-y))^(3/2)
    u[I,2] = 27*(x*(1-x))^(3/2)
  end
  u
end
xyd_brusselator = range(0,stop=1,length=32)

"""
2D Brusselator

```math
\\frac{\\partial u}{\\partial t} = 1 + u^2v - 4.4u + \\alpha(\\frac{\\partial^2 u}{\\partial x^2} + \\frac{\\partial^2 u}{\\partial y^2}) + f(x, y, t)
```
```math
\\frac{\\partial v}{\\partial t} = 3.4u - u^2v + \\alpha(\\frac{\\partial^2 u}{\\partial x^2} + \\frac{\\partial^2 u}{\\partial y^2})
```

where

``f(x, y, t) = \\begin{cases} 5 & \\quad \\text{if } (x-0.3)^2+(y-0.6)^2 ≤ 0.1^2 \\text{ and } t ≥ 1.1 \\\\ 0 & \\quad \\text{else} \\end{cases}``

and the initial conditions are

```math
u(x, y, 0) = 22\\cdot y(1-y)^{3/2}
```
```math
v(x, y, 0) = 27\\cdot x(1-x)^{3/2}
```

with the periodic boundary condition

```math
u(x+1,y,t) = u(x,y,t)
```
```math
u(x,y+1,t) = u(x,y,t)
```

From Hairer Norsett Wanner Solving Ordinary Differential Equations II - Stiff and Differential-Algebraic Problems Page 152
"""
prob_ode_brusselator_2d = ODEProblem(brusselator_2d_loop,
                                     init_brusselator_2d(xyd_brusselator),
                                     (0.,11.5),
                                     (3.4, 1., 10.,
                                      xyd_brusselator, step(xyd_brusselator),
                                      length(xyd_brusselator)))

const N_brusselator_1d = 40
const D_brusselator_u = CenteredDifference{Float64}(2,2,1/(N_brusselator_1d-1),
                                                    N_brusselator_1d)
const D_brusselator_v = CenteredDifference{Float64}(2,2,1/(N_brusselator_1d-1),
                                                    N_brusselator_1d)
function brusselator_1d(du, u_, p, t)
    A, B, α, buffer = p
    u = @view(u_[:, 1])
    v = @view(u_[:, 2])
    mul!(buffer, D_brusselator_u, u)
    Du = buffer
    @. du[:, 1] = A + u^2*v - (B+1)*u + α*Du

    mul!(buffer, D_brusselator_v, v)
    Dv = buffer
    @. du[:, 2] = B*u - u^2*v + α*Dv
    nothing
end
function init_brusselator_1d(N)
  u = zeros(N, 2)
  x = range(0, stop=1, length=N)
  for i in 1:N
    u[i, 1] = 1 + sin(2pi*x[i])
    u[i, 2] = 3.
  end
  u
end

"""
1D Brusselator

```math
\\frac{\\partial u}{\\partial t} = A + u^2v - (B+1)u + \\alpha\\frac{\\partial^2 u}{\\partial x^2}
```
```math
\\frac{\\partial v}{\\partial t} = Bu - u^2v + \\alpha\\frac{\\partial^2 u}{\\partial x^2}
```

and the initial conditions are

```math
u(x,0) = 1+\\sin(2π x)
```
```math
v(x,0) = 3
```

with the boundary condition

```math
u(0,t) = u(1,t) = 1
```
```math
v(0,t) = v(1,t) = 3
```

From Hairer Norsett Wanner Solving Ordinary Differential Equations II - Stiff and Differential-Algebraic Problems Page 6
"""
prob_ode_brusselator_1d = ODEProblem(brusselator_1d,
                                    init_brusselator_1d(N_brusselator_1d),
                                    (0.,10.),
                                    (1., 3., 1/50, zeros(N_brusselator_1d)))


text:
# ---
# title: 393. UTF-8 Validation
# id: problem393
# author: zhwang
# date: 2022-02-15
# difficulty: Medium
# categories: Bit Manipulation
# link: <https://leetcode.com/problems/utf-8-validation/description/>
# hidden: true
# ---
#
# A character in UTF8 can be from **1 to 4 bytes** long, subjected to the
# following rules:
#
#   1. For 1-byte character, the first bit is a 0, followed by its unicode code.
#   2. For n-bytes character, the first n-bits are all one's, the n+1 bit is 0, followed by n-1 bytes with most significant 2 bits being 10.
#
# This is how the UTF-8 encoding would work:
#
#
#
#        Char. number range  |        UTF-8 octet sequence
#           (hexadecimal)    |              (binary)
#        --------------------+---------------------------------------------
#        0000 0000-0000 007F | 0xxxxxxx
#        0000 0080-0000 07FF | 110xxxxx 10xxxxxx
#        0000 0800-0000 FFFF | 1110xxxx 10xxxxxx 10xxxxxx
#        0001 0000-0010 FFFF | 11110xxx 10xxxxxx 10xxxxxx 10xxxxxx
#
#
# Given an array of integers representing the data, return whether it is a valid
# utf-8 encoding.
#
# **Note:**
# The input is an array of integers. Only the **least significant 8 bits** of
# each integer is used to store the data. This means each integer represents
# only 1 byte of data.
#
# **Example 1:**
#
#
#
#     data = [197, 130, 1], which represents the octet sequence: **11000101 10000010 00000001**.
#
#     Return **true**.
#     It is a valid utf-8 encoding for a 2-bytes character followed by a 1-byte character.
#
#
# **Example 2:**
#
#
#
#     data = [235, 140, 4], which represented the octet sequence: **11101011 10001100 00000100**.
#
#     Return **false**.
#     The first 3 bits are all one's and the 4th bit is 0 means it is a 3-bytes character.
#     The next byte is a continuation byte which starts with 10 and that's correct.
#     But the second continuation byte does not start with 10, so it is invalid.
#
#
#
## @lc code=start
using LeetCode

## Automation
function valid_utf8(data::Vector{Int})::Bool
    ## 1 represents 0xxx.., 2 represents 10xx.., 3,4,5 similarly, 6 represents error
    switch(num::Int) = num >= 248 ? 6 : findfirst('0', bitstring(Int16(num))[end-7:end])
    ## 1 represents start, 2 represents 2-byte(need one more byte), 3, 4 similarly, 5 represents fail
    state = 1
    state_shift = Dict{Int,Vector{Int}}(1 => [1, 5, 2, 3, 4, 5])
    for i in 2:4
        state_shift[i] = [5, i-1, 5, 5, 5, 5]
    end
    for num in data
        println(state, num)
        state = state_shift[state][switch(num)] ## new state
        println(state)
        state == 5 && return false
    end
    state == 1
end
## @lc code=end

## @lc test=start
@testset "393.utf-8-validation.jl" begin
    @test valid_utf8([197, 130, 1]) == true
    @test valid_utf8([235, 140, 4]) == false
    @test valid_utf8([255]) == false
    @test valid_utf8([237]) == false
    @test valid_utf8([0]) == true
end
## @lc test=end

## @lc info=start
# link: [solution to 393](https://leetcode-cn.com/problems/utf-8-validation/solution/393-utf-8-bian-ma-yan-zheng-by-rex-rs-2ss8/)
## @lc info=end


text:
module PhysicsPrimitives2D

export PP2D
const PP2D = PhysicsPrimitives2D

import StaticArrays
const SA = StaticArrays
import LinearAlgebra
const LA = LinearAlgebra

include("utils.jl")
include("position_and_orientation.jl")
include("shapes.jl")
include("collision_detection.jl")
include("collision_manifold.jl")
include("abstract_body.jl")
include("physics.jl")

end


text:
using PyPlot, transD_GP, Random, LinearAlgebra

## model
ρbg        = 1000
zfixed     = [-1e5,   0,   (20:20:200)...]
rho        = [1e12,   ρbg, ρbg*ones(10)...]
ianom      = 5
rho[ianom] = 10.
δ          = 1e-4  # fd step size in log10

## Gauss Newton settings
use_fd          = false # in the one step for Gauss Newton
startfromcloser = false  # whether to start midway between anomaly and bg or plain bg in G-N step
λ²              = 5e-17 # model update damping

##  geometry and frequencies for VMD in frequency domain
nFreqsPerDecade     = 5
freqLowLimit        = 1e-3
freqHighLimit       = 1e5
freqs = 10 .^(log10(freqLowLimit):1/nFreqsPerDecade:log10(freqHighLimit))
rRx = 100.
zRx = -37.
zTx = -35.
nkᵣeval = 50
modelprimary = true
doconvramp = false
provideddt = true
Fvmd = transD_GP.AEM_VMD_HMD.HFieldDHT(
                      calcjacobian = true,
                      nmax = length(zfixed),
                      zTx    = zTx,
                      rRx    = rRx,
                      freqs  = freqs,
                      zRx    = zRx,
                      nkᵣeval = nkᵣeval,
                      modelprimary = modelprimary,
                      doconvramp = doconvramp,
                      provideddt = provideddt)
transD_GP.AEM_VMD_HMD.getfieldTD!(Fvmd, zfixed, rho)

## plot
times = Fvmd.interptimes
f, ax = plt.subplots(1,3, sharey=true, figsize=(8,6))
delz, delt = diff(zfixed[2:end]), diff(log10.(freqs))
tlabel = [(log10.(times) .- delt[1]/2)..., log10.(times)[end]+delt[end]/2]
zlab = [zfixed[2:end]..., zfixed[end]+delz[end]]
ax[1].step(log10.(rho[2:end]), zfixed[2:end])
im2 = ax[2].pcolormesh(tlabel, zlab, Fvmd.HTD_z_J_interp[2:end,:])
cb2 = colorbar(im2, ax=ax[2])
im3 = ax[3].pcolormesh(tlabel, zlab, log10.(abs.(Fvmd.HTD_z_J_interp[2:end,:])))
cb2 = colorbar(im3, ax=ax[3])
ax[1].set_xlabel("log10 Ω-m")
ax[1].invert_xaxis()
ax[2].invert_yaxis()
ax[2].set_xlabel("log10 s")
ax[3].set_xlabel("log10 s")
ax[1].set_ylabel("Depth m")
ax[2].set_title("∂f/∂(log10σ)")
ax[3].set_title("∂f/∂(log10σ)")

## approximate partials with central differences and compare
J = copy(Fvmd.HTD_z_J_interp)

function getJ_TD_noconv!(Fin, zin, ρ, Δ)
    # sensitivities computed in log10σ
    # same as in AEM_VMD_HMD code
    Jout = similar(Fin.HTD_z_J_interp)
    for i = 2:length(zin)
        modelfd = copy(ρ)
        modelfd[i] = 10 .^-(-log10(ρ[i]) .+ Δ/2)
        transD_GP.AEM_VMD_HMD.getfieldTD!(Fin, zin, modelfd)
        HTD_z1 = copy(Fin.HTD_z_interp)
        modelfd[i] = 10 .^-(-log10(ρ[i]) .- Δ/2)
        transD_GP.AEM_VMD_HMD.getfieldTD!(Fin, zin, modelfd)
        Jout[i,:] = (HTD_z1 - Fin.HTD_z_interp) / Δ
        end
    Jout
end
J_fd = getJ_TD_noconv!(Fvmd, zfixed, rho, δ)

## comparison plots
f, ax = plt.subplots(1,3, sharey=true, figsize=(10,6))
ax[1].step(log10.(rho[2:end]), zfixed[2:end])
im2 = ax[2].pcolormesh(tlabel, zlab, J[2:end,:])
cb2 = colorbar(im2, ax=ax[2])
im3 = ax[3].pcolormesh(tlabel, zlab, J_fd[2:end,:])
cb2 = colorbar(im3, ax=ax[3])
ax[1].set_xlabel("log10 Ω-m")
ax[1].invert_xaxis()
ax[2].invert_yaxis()
ax[2].set_xlabel("log10 Hz")
ax[3].set_xlabel("log10 s")
ax[1].set_ylabel("Depth m")
ax[2].set_title("∂f/∂(log10σ) analytic")
ax[3].set_title("∂f/∂(log10σ) FD")


text: retrieve(c::Client, ::Type{AuditLog}, g::AbstractGuild; kwargs...) = get_guild_audit_log(c, g.id; kwargs...)