text:
# test functions, potentially used by 2D and 3D tests

function filter_noise(A::AbstractArray)
  for i=1:length(A)
    if abs(A[i]) < 1e-8
      A[i] = 0
    end
  end

  return A
end

"""
  Test reverse mode of metric interpolation function
"""
function test_interp_rev(mesh)
  @testset "\n ----- Testing interp rev -----" begin

    dim = mesh.dim

    interp_data = PdePumiInterface.Interpolation(mesh)
    interp_datac = PdePumiInterface.Interpolation(mesh, Complex128)
    interp_data.bndry_arr[1] = Boundary(2, 1)
    interp_datac.bndry_arr[1] = Boundary(2, 1)

    if mesh.dim == 2
      dxidx_in = zeros(dim, dim, mesh.numNodesPerElement)
      jac_in = zeros(mesh.numNodesPerElement)

      for i=1:length(dxidx_in)
        dxidx_in[i] = i
      end
      for i=1:length(jac_in)
        jac_in[i] = i
      end
    else
      dxidx_in = mesh.dxidx[ :, :, :, 1]
      jac_in = mesh.jac[ :, 1]
    end

    dxidx_out = zeros(dim, dim, mesh.numNodesPerFace)
    jac_out = zeros(mesh.numNodesPerFace)

    PdePumiInterface.interpolateFace(interp_data, mesh.sbpface, dxidx_in, jac_in, dxidx_out, jac_out)

    # compute forward jacobian
    jac_size = dim*dim*mesh.numNodesPerFace
    jac_forward = zeros(jac_size, jac_size)

    dxidx_in2 = zeros(Complex128, size(dxidx_in)...)
    copy!(dxidx_in2, dxidx_in)
    jac_in2 = zeros(Complex128, size(jac_in)...)
    copy!(jac_in2, jac_in)

    dxidx_out2 = zeros(Complex128, dim, dim, mesh.numNodesPerFace)
    jac_out2 = zeros(Complex128, mesh.numNodesPerFace)

    h =1e-20
    pert = Complex128(0, h)

    for i=1:jac_size
      dxidx_in2[i] += pert
      PdePumiInterface.interpolateFace(interp_datac, mesh.sbpface, dxidx_in2, jac_in2, dxidx_out2, jac_out2)

      for j=1:jac_size
        jac_forward[j, i] = imag(dxidx_out2[j])/h
      end

      dxidx_in2[i] -= pert
    end

    # compute reverse mode
    jac_rev = zeros(jac_size, jac_size)
    dxidx_in_bar = zeros(dxidx_in)
    jac_in_bar = zeros(jac_in)
    dxidx_out_bar = zeros(dxidx_out)
    jac_out_bar = zeros(jac_out)

    for i=1:jac_size
      dxidx_out_bar[i] = 1
      fill!(dxidx_in_bar, 0.0)
      fill!(jac_in_bar, 0.0)
      fill!(jac_out_bar, 0.0)

      PdePumiInterface.interpolateFace_rev(interp_data, mesh.sbpface, dxidx_in,
                                       dxidx_in_bar, jac_in, jac_in_bar,
                                       dxidx_out, dxidx_out_bar, jac_out,
                                       jac_out_bar)

      for j=1:jac_size
        jac_rev[i, j] = dxidx_in_bar[j]
      end

      dxidx_out_bar[i] -= 1
    end

    @test isapprox( norm(jac_forward - jac_rev), 0.0) atol=1e-11

    # now test the full routine
    fill!(mesh.dxidx_bar, 0.0)
#    mesh.dxidx_face_bar[1] = 1
    fill!(mesh.dxidx_face_bar, 1.0)

    # check regular interfaces
    PdePumiInterface.interpolateMapping_rev(mesh)
    for i=1:mesh.numInterfaces
      el = mesh.interfaces[i].elementL
      for j=1:mesh.numNodesPerFace
        @test  norm(mesh.dxidx_bar[:, :, j, el])  > 0.0
      end
    end

    # check boundary faces
    fill!(mesh.dxidx_bar, 0.0)
    fill!(mesh.dxidx_bndry_bar, 1.0)
    PdePumiInterface.interpolateMapping_rev(mesh)
    for i=1:mesh.numBoundaryFaces
      el = mesh.bndryfaces[i].element
      for j=1:mesh.numNodesPerFace
        @test  norm(mesh.dxidx_bar[:, :, j, el])  > 0.0
      end
    end

    fill!(mesh.dxidx_bndry_bar, 1.0)

    # check sharedfaces
    for peer = 1:mesh.npeers
      fill!(mesh.dxidx_bar, 0.0)
      fill!(mesh.dxidx_sharedface[peer], 1.0)
      PdePumiInterface.interpolateMapping_rev(mesh)
      for i=1:mesh.numBoundaryFaces
        el = mesh.shared_interfacesfaces[p][i].elementL
        for j=1:mesh.numNodesPerFace
          @test  norm(mesh.dxidx_bar[:, :, j, el])  > 0.0
        end
      end

      fill!(mesh.dxidx_bar, 0.0)
      fill!(mesh.dxidx_sharedface[peer], 0.0)
    end

  end

end  # end function

"""
  Test the reverse mode of the curvilinear metric calculation.

  This function tests the reverse mode of calcEoneElement and then
  calls other functions to test the other parts of the calculation.
"""
function test_metrics_rev(mesh, mesh_c, sbp, opts)

  @testset "----- testing metric reverse mode -----" begin
    sbpface = mesh.sbpface
    nrm = rand(Complex128, mesh.dim, mesh.numNodesPerFace)
    nrm[:, :] = real(nrm)
    Rone = ones(mesh.numNodesPerFace)
    tmp = zeros(Complex128, mesh.numNodesPerFace)

    # test calcEoneElement

    Eone_el_c = zeros(Complex128, size(sbpface.perm, 1), mesh.dim)

    nin = length(nrm)
    nout = length(Eone_el_c)
    jac = zeros(nin, nout)
    jac2 = zeros(jac)

    h = 1e-20
    pert = Complex128(0, h)

    println("computing forward mode")
    for i=1:nin
      nrm[i] += pert
      fill!(Eone_el_c, 0.0)
      PdePumiInterface.calcEoneElement(sbpface, nrm, Rone, tmp, Eone_el_c)

      for j=1:nout
        jac[i, j] = imag(Eone_el_c[j])/h
      end

      nrm[i] -= pert
    end

    nrm_bar = zeros(mesh.dim, mesh.numNodesPerFace)
    tmp_bar = zeros(Float64, size(tmp)...)
    Eone_el_bar = zeros(Float64, size(Eone_el_c)...)

    println("computing reverse mode")
    for j=1:nout
      Eone_el_bar[j] = 1
      fill!(nrm_bar, 0.0)
      PdePumiInterface.calcEoneElement_rev(sbpface, nrm_bar, Rone, tmp_bar, Eone_el_bar)

      for i=1:nin
        jac2[i, j] = nrm_bar[i]
      end

      Eone_el_bar[j] = 0
    end

    @test isapprox( norm(jac - jac2)/length(jac2), 0.0) atol=1e-12

    # test assembleEone
    Eone_c = zeros(Complex128, mesh.numNodesPerElement, mesh.dim, 1)
    Eone_el_c[:, :] = real(Eone_el_c)
    elnum = 1
    facenum_local = 1

    nin = length(Eone_el_c)
    nout = length(Eone_c)

    jac = zeros(nin, nout)
    jac2 = zeros(jac)

    println("testing forward mode")
    for i=1:nin
      Eone_el_c[i] += pert
      fill!(Eone_c, 0.0)
      PdePumiInterface.assembleEone(sbpface, elnum, facenum_local, Eone_el_c, Eone_c)
      for j=1:nout
        jac[i, j] = imag(Eone_c[j])/h
      end

      Eone_el_c[i] -= pert
    end

    println("testing reverse mode")
    Eone_bar = zeros(Float64, mesh.numNodesPerElement, mesh.dim, 1)

    for j=1:nout
      Eone_bar[j] = 1

      fill!(Eone_el_bar, 0.0)
      PdePumiInterface.assembleEone_rev(sbpface, elnum, facenum_local, Eone_el_bar, Eone_bar)

      for i=1:nin
        jac2[i, j] = Eone_el_bar[i]
      end

      Eone_bar[j] = 0
    end

    @test isapprox( norm(jac - jac2), 0.0) atol=1e-12

    test_metric2_rev(mesh, mesh_c, sbp, opts)
    test_metrics3_rev(mesh, mesh_c, sbp, opts)
    test_metrics4_rev(mesh, mesh_c, sbp, opts)
    test_metrics5_rev(mesh, mesh_c, sbp, opts)

  end  # end facts block

  return nothing
end

"""
  Test calcEone_rev
"""
function test_metric2_rev(mesh, mesh_c, sbp, opts)

  @testset "----- Testing second part of Metric reverse mode -----" begin

    # test calcEone_rev
    # use finite differences because the perturbation is applied to the fields
    # of the mesh

    h = 1e-20
    pert = Complex128(0, h)

    nin = length(mesh.nrm_face) + length(mesh.nrm_bndry)
    for i=1:mesh.npeers
      nin += length(mesh.nrm_sharedface[i])
    end
    Eone = zeros(Complex128, mesh.numNodesPerElement, mesh.dim, mesh.numEl)
    nout = length(Eone)
#    nout = mesh.numNodesPerElement*mesh.dim*mesh.numEl

#    Eone_orig = zeros(Eone)

    jac = zeros(nin, nout)
    jac2 = zeros(jac)
    element_range = 1:mesh.numEl
 #   PdePumiInterface.calcEone(mesh, sbp, element_range, Eone_orig)

    println("forward mode")
    in_idx = 1
    for i=1:length(mesh_c.nrm_face)
      mesh_c.nrm_face[i] += pert
      fill!(Eone, 0.0)
      PdePumiInterface.calcEone(mesh_c, sbp, element_range, Eone)

      for j=1:nout
        jac[in_idx, j] = imag(Eone[j])/h  #(Eone[j] - Eone_orig[j])/pert
      end

      in_idx += 1
      mesh_c.nrm_face[i] -= pert
    end

    for i=1:length(mesh_c.nrm_bndry)
      mesh_c.nrm_bndry[i] += pert
      fill!(Eone, 0.0)
      PdePumiInterface.calcEone(mesh_c, sbp, element_range, Eone)
      for j=1:nout
        jac[in_idx, j] = imag(Eone[j])/h #(Eone[j] - Eone_orig[j])/pert
      end

      in_idx += 1
      mesh_c.nrm_bndry[i] -= pert
    end

    for peer=1:mesh_c.npeers
      nrm_peer = mesh_c.nrm_sharedface[peer]
      for i=1:length(nrm_peer)
        nrm_peer[i] += pert
        fill!(Eone, 0.0)
        PdePumiInterface.calcEone(mesh_c, sbp, element_range, Eone)
        for j=1:nout
          jac[in_idx, j] = imag(Eone[j])/h #(Eone[j] - Eone_orig[j])/pert
        end

        in_idx += 1
        nrm_peer[i] -= pert
      end
    end

    # reverse mode
    println("reverse mode")
    Eone_bar = zeros(mesh.numNodesPerElement, mesh.dim, mesh.numEl)
    for j=1:nout
      Eone_bar[j] = 1

      fill!(mesh.nrm_face_bar, 0.0)
      fill!(mesh.nrm_bndry_bar, 0.0)
      for i=1:mesh.npeers
        fill!(mesh.nrm_sharedface_bar[i], 0.0)
      end
      PdePumiInterface.calcEone_rev(mesh, sbp, element_range, Eone_bar)

      in_idx = 1
      for i=1:length(mesh.nrm_face_bar)
        jac2[in_idx, j] = mesh.nrm_face_bar[i]
        in_idx += 1
      end

      for i=1:length(mesh.nrm_bndry_bar)
        jac2[in_idx, j] = mesh.nrm_bndry_bar[i]
        in_idx += 1
      end

      for peer=1:mesh.npeers
        nrm_bar_peer = mesh.nrm_sharedface_bar[peer]
        for i=1:length(nrm_bar_peer)
          jac2[in_idx, j] = nrm_bar_peer[i]
          in_idx += 1
        end
      end

      Eone_bar[j] = 0

    end  # end loop j

    @test isapprox( maximum(abs.(jac - jac2)), 0.0) atol=1e-12
  end  # end facts block

  return nothing
end

"""
  Test getCurvilinearCoordinatesAndMetrics_rev (dxidx -> vert_coord, nrm)
"""
function test_metrics3_rev(mesh, mesh_c, sbp, opts)

  @testset "----- testing metrics reverse mode 3 -----" begin
#    nout = mesh.dim*mesh.numNodesPerElement*2
    nout = length(mesh.dxidx) + length(mesh.jac)
    nin = length(mesh.vert_coords) + length(mesh.nrm_face) + length(mesh.nrm_bndry)
    for i=1:mesh.npeers
      nin += length(mesh.nrm_sharedface[i])
    end
#    nin = length(mesh.vert_coords) + mesh.dim*mesh.numNodesPerElement*2

    jac = zeros(nin, nout)
    jac2 = zeros(jac)

    h = 1e-20
    pert = Complex128(0, h)

#    dxidx_orig = copy(mesh.dxidx)
#    jac_orig = copy(mesh.jac)
#    for i=1:mesh.npeers
#      nrm_sharedface[i] = copy(mesh.nrm_sharedface[i])
#    end

    # forward mode
    println("forward mode verts")
    in_idx = 1
    for i=1:length(mesh_c.vert_coords)
      mesh_c.vert_coords[i] += pert

      PdePumiInterface.getCurvilinearCoordinatesAndMetrics(mesh_c, sbp)

      out_idx = 1
      for j=1:length(mesh_c.dxidx)
        jac[in_idx, out_idx] = imag(mesh_c.dxidx[j])/h  #(mesh_c.dxidx[j] - dxidx_orig[j])/pert
        out_idx += 1
      end

      # TODO: uncomment when this is fixed
      for j=1:length(mesh_c.jac)
        #jac[in_idx, out_idx] = imag(mesh_c.jac[j])/h #(mesh_c.jac[j] - jac_orig[j])/pert
        out_idx += 1
      end

      in_idx += 1
      mesh_c.vert_coords[i] -= pert
    end

    for i=1:length(mesh_c.nrm_face)
      mesh_c.nrm_face[i] += pert
      PdePumiInterface.getCurvilinearCoordinatesAndMetrics(mesh_c, sbp)

      out_idx = 1
      for j=1:length(mesh_c.dxidx)
        jac[in_idx, out_idx] = imag(mesh_c.dxidx[j])/h #(mesh_c.dxidx[j] - dxidx_orig[j])/pert
        out_idx += 1
      end

      for j=1:length(mesh_c.jac)
        #jac[in_idx, out_idx] = imag(mesh_c.jac[j])/h  #(mesh_c.jac[j] - jac_orig[j])/pert
        out_idx += 1
      end

      in_idx += 1
      mesh_c.nrm_face[i] -= pert
    end

    for i=1:length(mesh_c.nrm_bndry)
      mesh_c.nrm_bndry[i] += pert
      PdePumiInterface.getCurvilinearCoordinatesAndMetrics(mesh_c, sbp)

      out_idx = 1
      for j=1:length(mesh_c.dxidx)
        jac[in_idx, out_idx] = imag(mesh_c.dxidx[j])/h  #(mesh_c.dxidx[j] - dxidx_orig[j])/pert
        out_idx += 1
      end

      # TODO: uncomment when this is fixed
      for j=1:length(mesh_c.jac)
        #jac[in_idx, out_idx] = imag(mesh_c.jac[j])/h #(mesh_c.jac[j] - jac_orig[j])/pert
        out_idx += 1
      end

      in_idx += 1
      mesh_c.nrm_bndry[i] -= pert
    end

    for peer=1:mesh_c.npeers
      println("peer ", peer)
      nrm_peer = mesh_c.nrm_sharedface[peer]
      for i=1:length(nrm_peer)
        nrm_peer[i] += pert
        PdePumiInterface.getCurvilinearCoordinatesAndMetrics(mesh_c, sbp)

        out_idx = 1
        for j=1:length(mesh_c.dxidx)
          jac[in_idx, out_idx] = imag(mesh_c.dxidx[j])/h #(mesh_c.dxidx[j] - dxidx_orig[j])/pert
          out_idx += 1
        end

        # uncomment when this is fixed
        for j=1:length(mesh_c.jac)
          #jac[in_idx, out_idx] = imag(mesh_c.jac[j])/h  #(mesh_c.jac[j] - jac_orig[j])/pert
          out_idx += 1
        end

        in_idx += 1
        nrm_peer[i] -= pert
      end
    end

    # reverse mode
    println("running reverse mode")
    fill!(mesh.dxidx_bar, 0.0)
    fill!(mesh.jac_bar, 0.0)
    out_idx = 1
    for j=1:length(mesh.dxidx)
      mesh.dxidx_bar[j] = 1
#      Eone_bar[j] = 1

      fill!(mesh.vert_coords_bar, 0.0)
      fill!(mesh.nrm_face_bar, 0.0)
      fill!(mesh.nrm_bndry_bar, 0.0)
      for i=1:mesh.npeers
        fill!(mesh.nrm_sharedface_bar[i], 0.0)
      end
      fill!(mesh.jac_bar, 0.0)
      PdePumiInterface.getCurvilinearCoordinatesAndMetrics_rev(mesh, sbp)

      in_idx = 1
      for i=1:length(mesh.vert_coords)
        jac2[in_idx, out_idx] = mesh.vert_coords_bar[i]
        in_idx += 1
      end

      for i=1:length(mesh.nrm_face)
        jac2[in_idx, out_idx] = mesh.nrm_face_bar[i]
        in_idx += 1
      end

      for i=1:length(mesh.nrm_bndry)
        jac2[in_idx, out_idx] = mesh.nrm_bndry_bar[i]
        in_idx += 1
      end

      for peer=1:mesh.npeers
        nrm_peer_bar = mesh.nrm_sharedface_bar[peer]
        for i=1:length(nrm_peer_bar)
          jac2[in_idx, out_idx] = nrm_peer_bar[i]
          in_idx += 1
        end
      end

      out_idx += 1
      mesh.dxidx_bar[j] = 0
    end  # end loop j

    fill!(mesh.jac_bar, 0.0)
    fill!(mesh.dxidx_bar, 0.0)
    for j=1:length(mesh.jac)
      #TODO: uncomment when this is fixed
      #mesh.jac_bar[j] = 1

      fill!(mesh.vert_coords_bar, 0.0)
      fill!(mesh.nrm_face_bar, 0.0)
      fill!(mesh.nrm_bndry_bar, 0.0)
      for i=1:mesh.npeers
        fill!(mesh.nrm_sharedface_bar[i], 0.0)
      end
      PdePumiInterface.getCurvilinearCoordinatesAndMetrics_rev(mesh, sbp)

      in_idx = 1
      for i=1:length(mesh.vert_coords)
        jac2[in_idx, out_idx] = mesh.vert_coords_bar[i]
        in_idx += 1
      end

      for i=1:length(mesh.nrm_face)
        jac2[in_idx, out_idx] = mesh.nrm_face_bar[i]
        in_idx += 1
      end

      for i=1:length(mesh.nrm_bndry)
        jac2[in_idx, out_idx] = mesh.nrm_bndry_bar[i]
        in_idx += 1
      end

      for peer=1:mesh.npeers
        nrm_peer_bar = mesh.nrm_sharedface_bar[peer]
        for i=1:length(nrm_peer_bar)
          jac2[in_idx, out_idx] = nrm_peer_bar[i]
          in_idx += 1
        end
      end

      out_idx += 1
      mesh.jac_bar[j] = 0
    end  # end loop j

    println("max diff = ", maximum(abs.(jac - jac2)))
    @test isapprox( maximum(abs.(jac - jac2)), 0.0) atol=1e-13
  end

  return nothing
end

function test_metrics4_rev(mesh, mesh_c, sbp, opts)

  # zero out all bar variables, just in case
  fill!(mesh.dxidx_bar, 0.0)
  fill!(mesh.jac_bar, 0.0)
  fill!(mesh.vert_coords_bar, 0.0)
  fill!(mesh.nrm_face_bar, 0.0)
  fill!(mesh.nrm_bndry_bar, 0.0)
  for i=1:mesh.npeers
    fill!(mesh.nrm_sharedface_bar[i], 0.0)
  end

  @testset "----- Testing metrics4_rev -----" begin
    nin = length(mesh.vert_coords)
    nout = length(mesh.nrm_bndry) + length(mesh.nrm_face) + length(mesh.coords_bndry)
    for i=1:mesh.npeers
      nout += length(mesh.nrm_sharedface[i])
    end

    jac = zeros(nin, nout)
    jac2 = zeros(jac)

    PdePumiInterface.getFaceCoordinatesAndNormals(mesh, sbp)
    nrm_bndry_orig = copy(mesh.nrm_bndry)
    nrm_face_orig = copy(mesh.nrm_face)
    nrm_sharedface_orig = copy(mesh.nrm_sharedface)
    for i=1:mesh.npeers
      nrm_sharedface_orig[i] = copy(mesh.nrm_sharedface[i])
    end

    h = 1e-20
    pert = Complex128(0, h)

    # forward mode
    in_idx = 1
    for i=1:length(mesh_c.vert_coords)
      mesh_c.vert_coords[i] += pert

      PdePumiInterface.getFaceCoordinatesAndNormals(mesh_c, sbp)

      out_idx = 1

      for j=1:length(mesh_c.nrm_bndry)
        jac[in_idx, out_idx] = imag(mesh_c.nrm_bndry[j])/h  #(mesh_c.nrm_bndry[j] - nrm_bndry_orig[j])/pert
        out_idx += 1
      end

      for j=1:length(mesh_c.nrm_face)
        jac[in_idx, out_idx] = imag(mesh_c.nrm_face[j])/h #(mesh_c.nrm_face[j] - nrm_face_orig[j])/pert
        out_idx += 1
      end

      for j=1:length(mesh_c.coords_bndry)
        jac[in_idx, out_idx] = imag(mesh_c.coords_bndry[j])/h
        out_idx += 1
      end

      for peer=1:mesh_c.npeers
        for j=1:length(mesh_c.nrm_sharedface[peer])
          jac[in_idx, out_idx] = imag(mesh_c.nrm_sharedface[peer][j])/h #(mesh_c.nrm_sharedface[peer][j] - nrm_sharedface_orig[peer][j])/pert
          out_idx += 1
        end
      end

      in_idx += 1
      mesh_c.vert_coords[i] -= pert
    end

    # reverse mode
    println("reverse mode")
    out_idx = 1
    for j=1:length(mesh.nrm_bndry)
      mesh.nrm_bndry_bar[j] = 1
      fill!(mesh.vert_coords_bar, 0.0)

      PdePumiInterface.getFaceCoordinatesAndNormals_rev(mesh, sbp)

      in_idx = 1
      for i=1:length(mesh.vert_coords_bar)
        jac2[in_idx, out_idx] = mesh.vert_coords_bar[i]
        in_idx += 1
      end

      out_idx += 1
      mesh.nrm_bndry_bar[j] = 0
    end

    for j=1:length(mesh.nrm_face)
      mesh.nrm_face_bar[j] = 1
      fill!(mesh.vert_coords_bar, 0.0)

      PdePumiInterface.getFaceCoordinatesAndNormals_rev(mesh, sbp)

      in_idx = 1
      for i=1:length(mesh.vert_coords_bar)
        jac2[in_idx, out_idx] = mesh.vert_coords_bar[i]
        in_idx += 1
      end

      out_idx += 1
      mesh.nrm_face_bar[j] = 0
    end

    for j=1:length(mesh.coords_bndry)
      mesh.coords_bndry_bar[j] = 1
      fill!(mesh.vert_coords_bar, 0.0)

      PdePumiInterface.getFaceCoordinatesAndNormals_rev(mesh, sbp)
      in_idx = 1
      for i=1:length(mesh.vert_coords_bar)
        jac2[in_idx, out_idx] = mesh.vert_coords_bar[i]
        in_idx += 1
      end

      out_idx += 1
      mesh.coords_bndry_bar[j] = 0
    end

    for peer=1:mesh.npeers
      nrm_sharedface_bar = mesh.nrm_sharedface_bar[peer]
      for j=1:length(nrm_sharedface_bar)
        nrm_sharedface_bar[j] = 1
        fill!(mesh.vert_coords_bar, 0.0)

        PdePumiInterface.getFaceCoordinatesAndNormals_rev(mesh, sbp)

        in_idx = 1
        for i=1:length(mesh.vert_coords_bar)
          jac2[in_idx, out_idx] = mesh.vert_coords_bar[i]
          in_idx += 1
        end

        out_idx += 1
        nrm_sharedface_bar[j] = 0
      end
    end

    diffnorm = maximum(abs.(jac - jac2))
    println("diffnorm = ", diffnorm)
    @test isapprox( diffnorm, 0.0) atol=1e-13
  end

  return nothing
end

"""
  Test everything together
"""
function test_metrics5_rev(mesh, mesh_c, sbp, opts)

  @testset "----- testing metrics5_rev -----" begin

    nin = length(mesh.vert_coords)

    nout = length(mesh.dxidx) + length(mesh.jac) + length(mesh.nrm_face) + length(mesh.nrm_bndry)
    for i=1:mesh.npeers
      nout += length(mesh.nrm_sharedface[i])
    end

    h = 1e-20
    pert = Complex128(0, h)

    jac = zeros(nin, nout)
    jac2 = zeros(jac)

    println("forward mode")
    for i=1:length(mesh_c.vert_coords)
      mesh_c.vert_coords[i] += pert
      PdePumiInterface.getFaceCoordinatesAndNormals(mesh_c, sbp)
      PdePumiInterface.getCurvilinearCoordinatesAndMetrics(mesh_c, sbp)

      out_idx = 1
      for j=1:length(mesh_c.dxidx)
        jac[i, out_idx] = imag(mesh_c.dxidx[j])/h
        out_idx += 1
      end

      for j=1:length(mesh_c.jac)
        jac[i, out_idx] = imag(mesh_c.jac[j])/h
        out_idx += 1
      end

      for j=1:length(mesh_c.nrm_face)
        jac[i, out_idx] = imag(mesh_c.nrm_face[j])/h
        out_idx += 1
      end

      for j=1:length(mesh_c.nrm_bndry)
        jac[i, out_idx] = imag(mesh_c.nrm_bndry[j])/h
        out_idx += 1
      end

      for peer=1:mesh_c.npeers
        nrm_peer = mesh_c.nrm_sharedface[peer]
        for j=1:length(nrm_peer)
          jac[i, out_idx] = imag(nrm_peer[j])/h
          out_idx += 1
        end
      end

      mesh_c.vert_coords[i] -= pert
    end

    println("reverse mode")
    zeroBarArrays(mesh)
    out_idx = 1
    for j=1:length(mesh.dxidx)
      mesh.dxidx_bar[j] = 1

      PdePumiInterface.getCurvilinearCoordinatesAndMetrics_rev(mesh, sbp)
      PdePumiInterface.getFaceCoordinatesAndNormals_rev(mesh, sbp)

      for i=1:length(mesh.vert_coords)
        jac2[i, out_idx] = mesh.vert_coords_bar[i]
      end

      out_idx += 1
      zeroBarArrays(mesh)
    end

    for j=1:length(mesh.jac)
      mesh.jac_bar[j] = 1  #TODO: uncomment when this is fixed

      #TODO: should these be in the other order?
#      getAllCoordinatesAndMetrics_rev(mesh, sbp)
      PdePumiInterface.getCurvilinearCoordinatesAndMetrics_rev(mesh, sbp)
      PdePumiInterface.getFaceCoordinatesAndNormals_rev(mesh, sbp)

      for i=1:length(mesh.vert_coords)
        jac2[i, out_idx] = mesh.vert_coords_bar[i]
      end

      out_idx += 1
      zeroBarArrays(mesh)
    end

    for j=1:length(mesh.nrm_face)
      mesh.nrm_face_bar[j] = 1

#      getAllCoordinatesAndMetrics_rev(mesh, sbp)
      PdePumiInterface.getCurvilinearCoordinatesAndMetrics_rev(mesh, sbp)
      PdePumiInterface.getFaceCoordinatesAndNormals_rev(mesh, sbp)

      for i=1:length(mesh.vert_coords)
        jac2[i, out_idx] = mesh.vert_coords_bar[i]
      end

      out_idx += 1
      zeroBarArrays(mesh)
    end

    for j=1:length(mesh.nrm_bndry)
      mesh.nrm_bndry_bar[j] = 1

#      getAllCoordinatesAndMetrics_rev(mesh, sbp)
      PdePumiInterface.getCurvilinearCoordinatesAndMetrics_rev(mesh, sbp)
      PdePumiInterface.getFaceCoordinatesAndNormals_rev(mesh, sbp)

      for i=1:length(mesh.vert_coords)
        jac2[i, out_idx] = mesh.vert_coords_bar[i]
      end

      out_idx += 1
      zeroBarArrays(mesh)
    end

    for peer=1:mesh.npeers
      nrm_bar_peer = mesh.nrm_sharedface_bar[peer]
      for j=1:length(nrm_bar_peer)
        nrm_bar_peer[j] = 1

#        getAllCoordinatesAndMetrics_rev(mesh, sbp)
        PdePumiInterface.getCurvilinearCoordinatesAndMetrics_rev(mesh, sbp)
        PdePumiInterface.getFaceCoordinatesAndNormals_rev(mesh, sbp)

        for i=1:length(mesh.vert_coords)
          jac2[i, out_idx] = mesh.vert_coords_bar[i]
        end

        out_idx += 1
        zeroBarArrays(mesh)
      end
    end

    @test isapprox( maximum(abs.(jac - jac2)), 0.0) atol=1e-12
  end

  return nothing
end

function test_coords_rev(mesh, sbp)

#  sbp_complex = TriSBP{Complex128}(degree=sbp.degree, reorder=sbp.reorder, internal=sbp.internal, vertices=sbp.vertices)

#  function coords_forward(vert_coords)
#    coords_it = vert_coords.'
#    coords_volume = SummationByParts.SymCubatures.calcnodes(sbp.cub, coords_it)
#    return coords_volume
# : end
  numVertsPerElement = mesh.numTypePerElement[1]
  # define indexing for input and output dimensions
  idx1(dim, vert, el) = dim + mesh.dim*(vert-1) + (el-1)*numVertsPerElement*mesh.dim
  idx2(dim, node, el) = dim + mesh.dim*(node-1) + (el-1)*mesh.numNodesPerElement*mesh.dim

  @testset "----- Testing coordiate reverse mode -----" begin
  # compute the jacobian with forward mode
  nin = mesh.dim*numVertsPerElement*mesh.numEl
  nout = mesh.dim*mesh.numNodesPerElement*mesh.numEl
  jac = zeros(nin, nout)

  vertcoords_in = zeros(Complex128, size(mesh.vert_coords)...)
#  coords_out = zeros(Complex128, size(mesh.coords)...)

  copy!(vertcoords_in, mesh.vert_coords)
  h = 1e-20
  pert = Complex128(0, h)
  for el=1:mesh.numEl
    for vert=1:numVertsPerElement
      for d=1:mesh.dim
        idx_in = idx1(d, vert, el)

        vertcoords_in[d, vert, el] += pert

        vertcoords_el = sview(vertcoords_in, :, :, el)
        coords_el = PdePumiInterface.vertToVolumeCoords(mesh, sbp, vertcoords_el)

        vertcoords_in[d, vert, el] -= pert

        # take advantage of the fact that the operation is element local
        for node_out=1:mesh.numNodesPerElement
          for d2=1:mesh.dim
            idx_out = idx1(d2, node_out, el)
            jac[idx_out, idx_in] = imag(coords_el[d2, node_out])/h
          end
        end

      end
    end
  end

  jac2 = zeros(jac)
  coords_bar = zeros(mesh.dim, mesh.numNodesPerElement, mesh.numEl)
  vertcoords_bar = zeros(size(mesh.vert_coords))
  # construct reverse mode jacobian
  for el=1:mesh.numEl
    for node=1:mesh.numNodesPerElement
      for d=1:mesh.dim
        coords_bar[d, node, el] = 1
        PdePumiInterface.volumeToVertCoords_rev(mesh, sbp, vertcoords_bar, coords_bar)

        coords_bar[d, node, el] = 0
        idx_out = idx1(d, node, el)
        for vert_in = 1:numVertsPerElement
          for d2=1:mesh.dim
            idx_in = idx1(d2, vert_in, el)
            jac2[idx_out, idx_in] = vertcoords_bar[d2, vert_in, el]
          end
        end

        fill!(vertcoords_bar, 0.0)

      end
    end
  end

  for i=1:nout
    for j=1:nin
      @test isapprox( jac[j, i], jac2[j, i]) atol=1e-13
    end
  end

  end

  return nothing
end

function test_submesh()

  @testset "----- Testing SubMesh -----" begin
    order = 1
    sbp = getTriSBPOmega(degree=order)
    vtx = sbp.vtx
    sbpface = TriFace{Float64}(order, sbp.cub, vtx)

    opts = PdePumiInterface.get_defaults()
    opts["order"] = order
    opts["coloring_distance"] = 2
    opts["dmg_name"] = ".null"
    opts["smb_name"] = "tri8l.smb"
    opts["numBC"] = 1
    opts["BC1"] = [0, 1, 2, 3]
    opts["BC1_name"] = "testBC"

    mesh = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=4)

    el_list = Cint[1, 2, 3, 4]

    submesh, subopts = PumiMeshDG2(mesh, sbp, opts, "resolveBC", el_list)

    @test ( submesh.order )== mesh.order
    @test ( submesh.numEl )== length(el_list)
    @test ( subopts["numBC"] )== 2
    @test ( subopts["BC2_name"] )== "resolveBC"
    @test ( (subopts["BC2"][1] in subopts["BC1"]) )== false

    # do a a dangling mesh
    el_list = Cint[1, 2, 5, 6]

    submesh, subopts = PumiMeshDG2(mesh, sbp, opts, "resolveBC", el_list)

    @test ( submesh.order )== mesh.order
    @test ( submesh.numEl )== length(el_list)
    @test ( subopts["numBC"] )== 2
    @test ( subopts["BC2_name"] )== "resolveBC"
    @test ( (subopts["BC2"][1] in subopts["BC1"]) )== false

    # test injection and rejection
    qold = zeros(Float64, mesh.numDof)
    qold2 = zeros(Float64, mesh.numDofPerNode, mesh.numNodesPerElement, mesh.numEl)
    for i=1:mesh.numEl
      for j=1:mesh.numNodesPerElement
        for k=1:mesh.numDofPerNode
          dof = mesh.dofs[k, j, i]
          qold[dof] = i + j + k
          qold2[k, j, i] = i + j + k
        end
      end
    end
    qnew = zeros(Float64, submesh.numDof)
    qnew2 = zeros(Float64, submesh.numDofPerNode, submesh.numNodesPerElement, submesh.numEl)

    injectionOperator(mesh, qold, submesh, qnew)
    injectionOperator(mesh, qold2, submesh, qnew2)

    for i=1:submesh.numDof
      qnew[i] += 1
      qnew2[i] += 1
    end

    qold3 = zeros(qold)
    qold4 = zeros(qold2)
    rejectionOperator(submesh, qnew, mesh, qold3)
    rejectionOperator(submesh, qnew2, mesh, qold4)

    # only the submesh elements got incremented
    for i=1:mesh.numEl
      for j=1:mesh.numNodesPerElement
        for k=1:mesh.numDofPerNode
          dof = mesh.dofs[k, j, i]
          if i in el_list
            @test ( qold3[dof] )== i + j + k + 1
            @test ( qold4[k, j, i] )== i + j + k + 1
          else
            @test ( qold3[dof] )== 0.0
            @test ( qold4[k, j, i] )== 0.0
          end
        end
      end
    end

    # test boundary interpolation array
    interp_arr = getBoundaryInterpArray(mesh, submesh)

    for iface in interp_arr
      @test ( iface.elementL in el_list )== true
      @test ( iface.elementR in el_list )== false
      @test ( iface.elementL == iface.elementR )== false
    end

  end  # end facts

end

function test_parallel_metrics(mesh, sbp, opts)

  # put known values into mesh arrays, do data exchange, check the values
  # came out correctly

  # need to skip any element that is on multiple boundaries
  skip_els = Set{Int}()
  for i=1:mesh.npeers
    el_list_i = mesh.local_element_lists[i]

    for j=1:(mesh.npeers-1)
      el_list_j = mesh.local_element_lists[i]
      common_els = intersect(el_list_i, el_list_j)
      for k=1:length(common_els)
        push!(skip_els, common_els[k])
      end
    end
  end
  # do one peer at a time in case a single element is shared with multiple peers
  @testset "----- Testing exchangeMetricInfo -----" begin
    for peer=1:mesh.npeers
      # vert_coords
      numel = mesh.local_element_counts[peer]
      rank = mesh.myrank

      for i=1:numel
        elnum = mesh.local_element_lists[peer]

        if elnum in skip_els
          skip_i = true
        else
          skip_i = false
        end

        val = 1 + rank
        for j=1:mesh.coord_numNodesPerElement
          for k=1:mesh.dim
            if skip_i
              tmp = -1
            else
              tmp = val
            end
            mesh.vert_coords[k, j, elnum] = tmp
            val += 1
          end
        end

        val = 2 + rank
        for j=1:mesh.numNodesPerElement
          for k=1:mesh.dim
            if skip_i
              tmp = -1
            else
              tmp = val
            end
            mesh.coords[k, j, elnum] = tmp
            val += 1
          end
        end

        val = 3 + rank
        for j=1:mesh.numNodesPerElement
          for d2=1:mesh.dim
            for d1=1:mesh.dim
              if skip_i
                tmp = -1
              else
                tmp = val
              end
              mesh.dxidx[d1, d2, j, elnum] = tmp
              val += 1
            end
          end
        end

        val = 4 + rank
        for j=1:mesh.numNodesPerElement
          if skip_i
            tmp = -1
          else
            tmp = val
          end
          mesh.jac[j, elnum] = tmp
          val += 1
        end
      end  # end loop i

    end  # end loop peer

    # do the data exchange
    PdePumiInterface.exchangeMetricInfo(mesh, sbp)

    # check that is came out right on the other end
    for peer=1:mesh.npeers
      numel = mesh.remote_element_counts[peer]
      obj = mesh.remote_metrics[peer]
      rank = mesh.peer_parts[peer]

      for i=1:numel

        val = 1 + rank
        for j=1:mesh.coord_numNodesPerElement
          for k=1:mesh.dim
            tmp = obj.vert_coords[k, j, i]
            if tmp >= 0
              @test ( tmp )== val
            end
            val += 1
          end
        end

        val = 2 + rank
        for j=1:mesh.numNodesPerElement
          for k=1:mesh.dim
            tmp = obj.coords[k, j, i]
            if tmp >= 0
              @test ( tmp )== val
            end
            val += 1
          end
        end

        val = 3 + rank
        for j=1:mesh.numNodesPerElement
          for d2=1:mesh.dim
            for d1=1:mesh.dim
              tmp = obj.dxidx[d1, d2, j, i]
              if tmp >= 0
                @test ( tmp )== val
              end
              val += 1
            end
          end
        end

        val = 4 + rank
        for j=1:mesh.numNodesPerElement
          tmp = obj.jac[j, i]
          if tmp >= 0
            @test ( tmp )== val
          end
          val += 1
        end

      end  # end loop i

    end  # end loop peer

  end  # end facts block

  return nothing
end

function testSurfaceNumbering(mesh, sbp, opts)

  @testset "----- Testing Surface Numbering -----" begin
    nBCs = opts["numBC"]
    bc_nums =  [nBCs]  # only do the last BC  #TODO: re-enable this
#    bc_nums = 1:nBCs
    numFacePts, n_face, face_verts = numberSurfacePoints(mesh, bc_nums)

    @test ( length(face_verts) )== numFacePts

    geo_tags = Array{Int}(0)
    for bc in bc_nums
      geo_tags_bc = opts["BC$bc"]
      for j=1:length(geo_tags_bc)
        push!(geo_tags, geo_tags_bc[j])
      end
    end

    # check that if the MeshEntity is classified in an edge/face, then it
    # is part of the geometry
    # this is a bit too loose because some MeshEntities are classified on
    # the geometric entities that are on the closure of the specified geometry
    for vert in face_verts
      me =  apf.toModel(mesh.m_ptr, vert)
      me_type = apf.getModelType(mesh.m_ptr, me)  # dimension of model entity
      me_tag = apf.getModelTag(mesh.m_ptr, me)

      if me_type == mesh.dim - 1
        @test ( me_tag in geo_tags )== true
      end

    end

    # check that all entities with number numFacePts + 1 are not on the geometry

    seen_nums = zeros(Int, numFacePts)
    entities = Array{Vector{Ptr{Void}}}(2)
    entities[1] = mesh.verts
    entities[2] = mesh.edges
    for edim = 1:length(entities)
      if apf.hasNodesIn(mesh.coordshape_ptr, edim - 1)
        for vert in entities[edim]
          n_v = apf.getNumberJ(n_face, vert, 0, 0)
          me =  apf.toModel(mesh.m_ptr, vert)
          me_type = apf.getModelType(mesh.m_ptr, me)
          me_tag = apf.getModelTag(mesh.m_ptr, me)

          if n_v > numFacePts
            if me_type == mesh.dim - 1
              @test ( !(me_tag in geo_tags) )== true
            end
          else
            seen_nums[n_v] += 1  # for uniqueness check
            if me_type == mesh.dim - 1
              @test ( me_tag in geo_tags )== true
            end
          end
        end  # end loop vert
      end  # end if
    end  # end loop edim

    @test ( maximum(seen_nums) )== 1
    @test ( minimum(seen_nums) )== 1
  end  # end facts block

  return nothing
end

"""
  Do sanity checks on the coordinate field number.  Currently tests

   * all coordinate field dofs have numbers assigned
   * all the numbers are in the range 1 to mesh.dim * coord_numNodesPerElement
   * all numbers from 1 to coord_numNodesPerElement exist exactly once
"""
function test_coordNumbering(mesh)

  @testset "Testing coordinate field Numbering" begin

    nums = zeros(mesh.dim*mesh.coord_numNodes)
    # check verts
    for entity in mesh.verts
      for i=1:mesh.dim
        val = apf.getNumberJ(mesh.coord_nodenums_Nptr, entity, 0, i-1)
        @test val >= 1
        @test val <= mesh.dim*mesh.coord_numNodes

        nums[val] += 1
      end
    end

    if mesh.coord_order == 2
      for entity in mesh.edges
        for i=1:mesh.dim
          val = apf.getNumberJ(mesh.coord_nodenums_Nptr, entity, 0, i-1)
          @test val >= 1
          @test val <= mesh.dim*mesh.coord_numNodes

          nums[val] += 1
        end
      end
    end

    # every number exists once
    for val in nums
      @test val == 1
    end

  end  # end testset

  return nothing
end

"""
  Test coordinate field transformation functions
"""
function test_coord_field(mesh::PumiMeshDG{T}) where {T}

  @testset "Testing coordinate field transformation" begin

    # test round tripping the vert_coords array
    coord_arr = copy(mesh.vert_coords)
    coord_vec = zeros(T, mesh.dim*mesh.coord_numNodes)

    coords3DTo1D(mesh, coord_arr, coord_vec, PdePumiInterface.AssignReduction{T}())
    fill!(coord_arr, 0.0)
    coords1DTo3D(mesh, coord_vec, coord_arr, PdePumiInterface.AssignReduction{T}())

    @test vecnorm(coord_arr - mesh.vert_coords) < 1e-13

  end

  return nothing
end

"""
  Returns an array of the specified size with random values for the real part
  and zeros for the imaginary part
"""
function rand_realpart(dims...)

  a = rand(Complex128, dims...)
  for i=1:length(a)
    a[i] = real(a[i])
  end

  return a
end

function rand_realpart!(A::AbstractArray)

  rand!(A)
  for i=1:length(A)
    A[i] = real(A[i])
  end

  return nothing
end

import PdePumiInterface: RemoteMetrics

"""
  Tests back propigating the metrics to the 1D vector form
"""
function test_metrics_rev_1d(mesh::PumiMesh{T}, sbp, opts) where {T}

  println("testing metrics_rev_1d")
  srand(1234)
  h = 1e-20
  pert = Complex128(0, h)

  xvec = zeros(Complex128, mesh.dim*mesh.coord_numNodes)
  xvec_dot = rand_realpart(size(xvec))
  xvec_bar = zeros(Complex128, length(xvec))

  dxidx_bar       = rand_realpart(size(mesh.dxidx))
  jac_bar         = rand_realpart(size(mesh.jac))
  nrm_bndry_bar   = rand_realpart(size(mesh.nrm_bndry))
  nrm_face_bar    = rand_realpart(size(mesh.nrm_face_bar))
  coords_bndry_bar    = rand_realpart(size(mesh.coords_bndry_bar))
  nrm_sharedface_bar = Array{Array{Complex128, 3}}(mesh.npeers)
  remote_metrics_bar = Array{RemoteMetrics{Complex128}}(mesh.npeers)
  for i=1:mesh.npeers
    nrm_sharedface_bar[i] = rand_realpart(size(mesh.nrm_sharedface[i]))
    remote_metrics_bar[i] = RemoteMetrics(mesh, i, islocal=false)
    obj = remote_metrics_bar[i]
    rand_realpart!(obj.vert_coords)
    #rand_realpart!(obj.coords)
    rand_realpart!(obj.dxidx)
    rand_realpart!(obj.jac)
  end

  zeroBarArrays(mesh)
  vert_coords_orig = copy(mesh.vert_coords)

  # get unique-ified xvec_dot
  coords3DTo1D(mesh, mesh.vert_coords, xvec, PdePumiInterface.AssignReduction{T}(), parallel=false)
  xvec .+= pert*xvec_dot
  coords1DTo3D(mesh, xvec, mesh.vert_coords, parallel=true)
  xvec_dot = imag(xvec)/h  # unique-ified version

  diff = maximum(abs.(vert_coords_orig - real(mesh.vert_coords)))
  @assert diff < 1e-15

  # test that coords3DTo1D is the reverse mode of coords1DTo3D
  vert_coords2 = zeros(mesh.vert_coords)
  vert_coords_bar = rand_realpart(size(mesh.vert_coords))
  xvec_bar = zeros(Complex128, length(xvec))

  # forward mode
  coords1DTo3D(mesh, xvec, vert_coords2, parallel=false)
  #xvec .-= pert*xvec_dot
  for i=1:length(xvec)
    xvec[i] = real(xvec[i])
  end
  val1 = sum(imag(vert_coords2)/h .* vert_coords_bar)

  # reverse mode
  coords3DTo1D(mesh, vert_coords_bar, xvec_bar, parallel=true)
  val2 = sum(xvec_bar .* xvec_dot)

  val1 = MPI.Allreduce(val1, MPI.SUM, mesh.comm)
  val2 = MPI.Allreduce(val2, MPI.SUM, mesh.comm)
  @test abs(val1 - val2) < max(abs(val1)*1e-13, 1e-12)

  #----------------------------------------------------------------------------
  # test back-propigation of the metrics (locally)
  zeroBarArrays(mesh)
  fill!(xvec_bar, 0)
  vert_coords_dot = rand_realpart(size(mesh.vert_coords))

  for i=1:length(mesh.vert_coords)
    mesh.vert_coords[i] = real(mesh.vert_coords[i])
  end

  mesh.vert_coords .+= pert*vert_coords_dot
  PdePumiInterface.recalcCoordinatesAndMetrics(mesh, sbp, opts)
  mesh.vert_coords .-= pert*vert_coords_dot

  val1 = sum(imag(mesh.dxidx)/h .* dxidx_bar)               +
         sum(imag(mesh.jac)/h .* jac_bar)                  +
         sum(imag(mesh.nrm_bndry)/h .* nrm_bndry_bar)      +
         sum(imag(mesh.nrm_face)/h .* nrm_face_bar)          +
         sum(imag(mesh.coords_bndry)/h .* coords_bndry_bar)
  for i=1:mesh.npeers
    val1 += sum(imag(mesh.nrm_sharedface[i])/h .* nrm_sharedface_bar[i])
  end

  for i=1:length(mesh.vert_coords)
    mesh.vert_coords[i] = real(mesh.vert_coords[i])
  end

  copy!(mesh.dxidx_bar, dxidx_bar)
  copy!(mesh.jac_bar, jac_bar)
  copy!(mesh.nrm_bndry_bar, nrm_bndry_bar)
  copy!(mesh.nrm_face_bar, nrm_face_bar)
  copy!(mesh.coords_bndry_bar, coords_bndry_bar)
  for i=1:mesh.npeers
    copy!(mesh.nrm_sharedface_bar[i], nrm_sharedface_bar[i])
  end

  getAllCoordinatesAndMetrics_rev(mesh, sbp, opts)
#  getAllCoordinatesAndMetrics_rev(mesh, sbp, opts, xvec_bar, parallel=true)
  val2 = sum(mesh.vert_coords_bar .* vert_coords_dot)
  #val2 = sum(xvec_bar .* xvec_dot)

  @test abs(val1 - val2) < max(abs(val1)*1e-12, 1e-12)

  #----------------------------------------------------------------------------
  # test back-propigation of metrics (parallel)

  zeroBarArrays(mesh)
  fill!(xvec_bar, 0)

  # forward mode
  xvec .+= pert*xvec_dot
  coords1DTo3D(mesh, xvec, mesh.vert_coords, parallel=true)
  PdePumiInterface.recalcCoordinatesAndMetrics(mesh, sbp, opts)
  xvec .-= pert*xvec_dot

  val1 = sum(imag(mesh.dxidx)/h .* dxidx_bar)               +
         sum(imag(mesh.jac)/h .* jac_bar)                  +
         sum(imag(mesh.nrm_bndry)/h .* nrm_bndry_bar)      +
         sum(imag(mesh.nrm_face)/h .* nrm_face_bar)          +
         sum(imag(mesh.coords_bndry)/h .* coords_bndry_bar)
  for i=1:mesh.npeers
    val1 += sum(imag(mesh.nrm_sharedface[i])/h .* nrm_sharedface_bar[i])

    arr_dot = mesh.remote_metrics[i].vert_coords
    arr_bar = remote_metrics_bar[i].vert_coords
    val1 += sum(imag(arr_dot)/h .* arr_bar)

    #arr_dot = mesh.remote_metrics[i].coords
    #arr_bar = remote_metrics_bar[i].coords
    #val1 += sum(imag(arr_dot)/h .* arr_bar)

    arr_dot = mesh.remote_metrics[i].dxidx
    arr_bar = remote_metrics_bar[i].dxidx
    val1 += sum(imag(arr_dot)/h .* arr_bar)

    arr_dot = mesh.remote_metrics[i].jac
    arr_bar = remote_metrics_bar[i].jac
    val1 += sum(imag(arr_dot)/h .* arr_bar)

  end

  # remove complex parts
  for i=1:length(mesh.vert_coords)
    mesh.vert_coords[i] = real(mesh.vert_coords[i])
  end
  PdePumiInterface.recalcCoordinatesAndMetrics(mesh, sbp, opts)

  # reverse mode
  copy!(mesh.dxidx_bar, dxidx_bar)
  copy!(mesh.jac_bar, jac_bar)
  copy!(mesh.nrm_bndry_bar, nrm_bndry_bar)
  copy!(mesh.nrm_face_bar, nrm_face_bar)
  copy!(mesh.coords_bndry_bar, coords_bndry_bar)
  for i=1:mesh.npeers
    copy!(mesh.nrm_sharedface_bar[i], nrm_sharedface_bar[i])
    copy!(mesh.remote_metrics_bar[i], remote_metrics_bar[i])
  end

  getAllCoordinatesAndMetrics_rev(mesh, sbp, opts, xvec_bar)
  val2 = sum(xvec_bar .* xvec_dot)
  val1 = MPI.Allreduce(val1, MPI.SUM, mesh.comm)
  val2 = MPI.Allreduce(val2, MPI.SUM, mesh.comm)

  @test abs(val1 - val2) < max(abs(val1)*1e-12, 1e-12)

  return nothing
end

function test_geoNums(mesh)

  @testset "testing geometric numbering" begin
    geonums = mesh.geoNums

    for dim=0:mesh.dim
      for e in apf.MeshIterator(mesh.m_ptr, dim)
        me = apf.toModel(mesh.m_ptr, e)
        me_dim = apf.getModelType(mesh.m_ptr, me)

        for j=1:mesh.coord_numNodesPerType[dim+1]

          for k=0:(me_dim-1)
            @assert apf.isNumbered(geonums.xiNums, e, j-1, k)
            n_k = apf.getNumberJ(geonums.xiNums, e, j-1, k)
            @test n_k <= geonums.numXiDofs
          end

          for k=me_dim:(mesh.dim-1)
            @assert apf.isNumbered(geonums.xiNums, e, j-1, k)
            n_k = apf.getNumberJ(geonums.xiNums, e, j-1, k)
            @test n_k == geonums.numXiDofs + 1
          end
        end   # end j
      end # end iterator
    end  # end dim

  end # end testset

  return nothing
end

function test_geoMapping(mesh)
# mesh must have loaded a CAD-based geometric model

  @testset "testing geoMapping" begin
    println("testing geometric mapping")
    xivec = zeros(mesh.geoNums.numXiDofs)
    coords_arr = copy(mesh.vert_coords)

    coords_xyzToXi(mesh, coords_arr, xivec)
    fill!(coords_arr, 0)
    coords_XiToXYZ(mesh, xivec, coords_arr)

    println("maxdiff = ", maximum(abs.(coords_arr - mesh.vert_coords)))
    @test maximum(abs.(coords_arr - mesh.vert_coords)) < 1e-8 # CAD is not that accurate

  end

  return nothing
end

function test_geoWrapping(mesh)
# mesh must have loaded a CAD-based geometric model

  println("testing geometric wrapping")
  @testset "testing geoWrapping" begin
    # test the kernel function
    rng = [0.0, 1.0]
    @test abs(PdePumiInterface.wrapPeriodicXi(rng, 1.5) - 0.5) < 1e-12
    @test abs(PdePumiInterface.wrapPeriodicXi(rng, 1.4) - 0.4) < 1e-12
    @test abs(PdePumiInterface.wrapPeriodicXi(rng, 2.4) - 0.4) < 1e-12
    @test abs(PdePumiInterface.wrapPeriodicXi(rng, 3.4) - 0.4) < 1e-12

    rng = [1.5, 3.5]
    @test abs(PdePumiInterface.wrapPeriodicXi(rng, 4.0) - 2.0) < 1e-12

    val = Complex128(4.0, 1.5)
    val2 = Complex128(2.0, 1.5)
    @test abs(PdePumiInterface.wrapPeriodicXi(rng, val) - val2) < 1e-12

    # test the overall function

    # get periodic model entities
    g = apf.getModel(mesh.m_ptr)
    periodic_mes = Array{gmi.ModelEntity}(0)
    for me in gmi.Gmi_iter(g, mesh.dim-1)
      if gmi.periodic(g, me, 0)  # only check first dimension
        push!(periodic_mes, me)
      end
    end

    xvec = zeros(Float64, mesh.geoNums.numCoordDofs)
    op = PdePumiInterface.AssignReduction{Float64}()
    coords3DTo1D(mesh, mesh.vert_coords, xvec, op, parallel=false)

    xivec = zeros(mesh.geoNums.numXiDofs)
    coords_xyzToXi(mesh, xvec, xivec)
    xivec2 = copy(xivec)

    rng = zeros(Float64, 2)
    geonums = mesh.geoNums
    for (e, dim) in apf.FieldEntityIt(mesh.m_ptr, mesh.coordshape_ptr)
      me = apf.toModel(mesh.m_ptr, e)
      me_dim = apf.getModelType(mesh.m_ptr, me)
      if me in periodic_mes
        # add a multiple of the period to the xivec2
        gmi.range(g, me, 0, rng)
        delta = rng[2] - rng[1]

        idx = apf.getNumberJ(geonums.xiNums, e, 0, 0)
        xivec2[idx] += 2*delta
      end
    end

    # call wrapXiVals, check xivec2 == xivec
    PdePumiInterface.wrapXiCoords(mesh, xivec2)
    @test maximum(abs.(xivec - xivec2)) < 1e-12

  end  # end testset

  return nothing
end

function test_geoDerivative(mesh)
# mesh must have loaded a CAD-based geometric model.
# test against finite difference

  function f_x(i, xvec::AbstractVector)
    # compute f(x)

      coeff1 = i % 10
      coeff2 = i % 5

      x_i = xvec[i]
#      val = coeff1*x_i*x_i + coeff2*x_i + i
      val = coeff1*x_i

    return val
  end  # end function

  function f_x_deriv(i, xvec::AbstractVector, df_dx::AbstractVector)
    # compute f(x)

    #for i=1:length(xvec)
#    i = 19
      coeff1 = i % 10
      coeff2 = i % 5

      x_i = xvec[i]
      #df_dx[i] = 2*coeff1*x_i + coeff2
      df_dx[i] = coeff1
    #end
  end  # end function

  function f2_x(xvec::AbstractVector)

    val = 0.0
    for i=1:length(xvec)
      val += xvec[i]*xvec[i]
    end

    return val
  end

  function f2_x_deriv(xvec::AbstractVector, df_dx::AbstractVector)
    for i=1:length(xvec)
      df_dx[i] = 2*xvec[i]
    end
  end

  @testset "Geometric Derivative" begin
    println("\ntesting geometric derivative")
    xvec = getXCoords(mesh)
    xvec_pert = copy(xvec)
    xivec = getXiCoords(mesh)
    #=
    xivec = zeros(mesh.geoNums.numXiDofs)
    coords_xyzToXi(mesh, xvec, xivec)
    coords_XiToXYZ(mesh, xivec, xvec)  # hopefully this reduces numerical
                                       # error when compared to coords_XiToXYZ
    =#
    df_dx = zeros(mesh.geoNums.numCoordDofs)
    df_dxi = zeros(mesh.geoNums.numXiDofs)

    xi_indices = constructGeoMapping(mesh)
    h = 1e-6

#=
    # test each component individually
    # this test is rather slow, so don't run it.
    maxdiff = 0.0
    for i=1:length(xvec)
      # function at initial x value
      val1 = f_x(i, xvec)
      fill!(df_dx, 0)
      f_x_deriv(i, xvec, df_dx)
      coords_dXTodXi(mesh, df_dx, df_dxi)

      # only do j values related to i
      for j in xi_indices[i]
        if length(xi_indices[i]) != 1
          continue
        end
        # compute finite difference
        xivec[j] += h
        coords_XiToXYZ(mesh, xivec, xvec_pert)
        val2 = f_x(i, xvec_pert)
        xivec[j] -= h

        dfdxi_fd = (val2 - val1)/h
        #println("df_dxi = ", df_dxi[j])
        #println("df_dxi_fd = ", dfdxi_fd)

        diff = abs(df_dxi[j] - dfdxi_fd)
        #@test abs(df_dxi[j] - dfdxi_fd) < 1e-5

        if diff > maxdiff
          maxdiff = diff
        end
      end  # end j
    end  # end i
    println("maxdiff = ", maxdiff)
=#
    h = 1e-6
    pert = rand(length(xivec))

    # finite difference
    J1 = f2_x(xvec)
    for i=1:length(pert)
      xivec[i] += h*pert[i]
    end
    coords_XiToXYZ(mesh, xivec, xvec_pert)
    J2 = f2_x(xvec_pert)
    for i=1:length(pert)
      xivec[i] -= h*pert[i]
    end
    val1 = (J2 - J1)/h

    # using CAD derivative
    dJdx = zeros(xvec)
    dJdxi = zeros(xivec)
    f2_x_deriv(xvec, dJdx)
    coords_dXTodXi(mesh, dJdx, dJdxi)

    val2 = dot(dJdxi, pert)

    @test abs(val1 - val2) < 1e-2

  end  # end testset

  return nothing
end

function constructGeoMapping(mesh)
# get array of all xi indices that correspond to a given x index

  geonums = mesh.geoNums
  indices = Array{Array{Int, 1}}(geonums.numCoordDofs)
  for i=1:mesh.geoNums.numCoordDofs
    indices[i] = Array{Int}(0)
  end

  for (e, dim) in apf.FieldEntityIt(mesh.m_ptr, mesh.coordshape_ptr)
    for j=1:mesh.coord_numNodesPerType[dim+1]
      for k=1:mesh.dim
        idx_in = apf.getNumberJ(geonums.coordNums, e, j-1, k-1)
        for d=1:mesh.dim
          idx_out = apf.getNumberJ(geonums.xiNums, e, j-1, d-1)
          if idx_out != geonums.numXiDofs + 1
            push!(indices[idx_in], idx_out)
          end  # end if
        end  # end d
      end   # end k
    end  # end j
  end  # end iterator

  return indices
end

function test_setPoint(mesh, opts)

  coords_orig = zeros(Float64, 3)
  coords = zeros(Float64, 3)
  coords2 = zeros(Float64, 3)
  coords3 = zeros(Float64, 3)
  xi_orig = zeros(Float64, 3)
  xi = zeros(Float64, 3)
  xi3 = zeros(Float64, 3)
  rng = zeros(Float64, 2)
  geoname = opts["dmg_name"]
  if endswith(geoname, ".x_t") || endswith(geoname, ".smd")
    println("\nTesting setPoint")
    @test mesh.geoNums.can_eval

    println("testing setCoordsXi")
    for (e, dim) in apf.FieldEntityIt(mesh.m_ptr, mesh.coordshape_ptr)
      for j=1:mesh.coord_numNodesPerType[dim+1]
        g = apf.getModel(mesh.m_ptr)
        me = apf.toModel(mesh.m_ptr, e)
        me_dim = apf.getModelType(mesh.m_ptr, me)

        # skip things that don't have parametric coords/can't change them
        if me_dim == 0 || me_dim == 3
          continue
        end

        getCoordsXi(mesh, e, j-1, xi_orig, coords_orig)
        getCoordsXi(mesh, e, j-1, xi, coords)

        for i=1:me_dim
          gmi.range(g, me, i-1, rng)
          offset = 0.01
          if xi[i] < rng[2] - offset  # avoid out of range values
            xi[i] += 0.01
          end
        end

        setCoordsXi(mesh, e, j-1, xi, coords2)
        getCoordsXi(mesh, e, j-1, xi3, coords3)

        @test maximum(abs.(coords2 - coords3)) < 1e-8
        @test maximum(abs.(xi - xi3)) < 1e-13
      end
    end

  end  # end if

  # test setCoords
  println("testing setCoords")
  for (e, dim) in apf.FieldEntityIt(mesh.m_ptr, mesh.coordshape_ptr)
    for j=1:mesh.coord_numNodesPerType[dim+1]
      g = apf.getModel(mesh.m_ptr)
      me = apf.toModel(mesh.m_ptr, e)
      me_dim = apf.getModelType(mesh.m_ptr, me)

      getCoords(mesh, e, j-1, coords_orig)

      if mesh.geoNums.can_eval
        if me_dim > 0 && me_dim < 3
          getCoordsXi(mesh, e, j-1, xi_orig, coords)
        else
          getCoords(mesh, e, j-1, coords)
        end

        # test snapping
        for i=1:mesh.dim
          coords[i] += 0.1
        end

        setCoords(mesh, e, j-1, coords, true)

        if me_dim == 0
          getCoords(mesh, e, j-1, coords3)
          @test maximum(abs.(coords3 - coords_orig)) < 1e-8
        elseif me_dim == 3
          getCoords(mesh, e, j-1, coords3)
          @test maximum(abs.(coords3 - coords)) < 1e-8
        else
          getCoordsXi(mesh, e, j-1, xi3, coords3)
          gmi.geval(g, me, xi3, coords2)
          @test maximum(abs.(coords2 - coords)) < 1e-8
          @test maximum(abs.(coords3 - coords)) < 1e-8
        end

        # test non snapping
        fill!(coords2, 0); fill!(coords3, 0)
        fill!(xi, 0); fill!(xi3, 0)

        setCoords(mesh, e, j-1, coords_orig, true)
        for i=1:3
          coords[i] = coords_orig[i] + 0.1
        end
        setCoords(mesh, e, j-1, coords, false)

        if me_dim == 0
          getCoords(mesh, e, j-1, coords3)
          @test maximum(abs.(coords3 - coords)) < 1e-8
        elseif me_dim == 3
          getCoords(mesh, e, j-1, coords3)
          @test maximum(abs.(coords3 - coords)) < 1e-8
        else
          # check that the xi coordinates changed, even though they are
          # not consistent with the xyz coordinates
          getCoordsXi(mesh, e, j-1, xi3, coords3)
          gmi.geval(g, me, xi3, coords2)
          coords_tmp = zeros(Float64, 3)
          gmi.geval(g, me, xi_orig, coords_tmp)
          @test maximum(abs.(coords2 - coords)) > 1e-8
          @test maximum(abs.(coords3 - coords)) < 1e-8
        end
      else  # cannot eval

        for i=1:mesh.dim
          coords[i] += 0.1
        end

        setCoords(mesh, e, j-1, coords, true)
        getCoords(mesh, e, j-1, coords2)

        @test maximum(abs.(coords2 - coords)) < 1e-13
      end
    end
  end

  return nothing
end

"""
  Tests the single-element version of update_coords with snapping
"""
function test_geoWarp(mesh, sbp, opts)

  # test single-element update_coords
  coords_orig = copy(mesh.vert_coords)

  # double all coordinates
  coords_i = zeros(Float64, 2, 3)
  for i=1:mesh.numEl
    for j=1:3
      coords_i[1, j] = coords_orig[1, j, i] + 0.0001
      coords_i[2, j] = coords_orig[2, j, i] + 0.0001
    end
    update_coords(mesh, i, coords_i)
  end

  commit_coords(mesh, sbp, opts)

  # check xyz and xi are consistent
  xyz_j = zeros(Float64, 3)
  xi_j = zeros(Float64, 2)
  xyz2_j = zeros(Float64, 3)
  g = apf.getModel(mesh.m_ptr)
  for e in apf.MeshIterator(mesh.m_ptr, 0)
    apf.getPoint(mesh.m_ptr, e, 0, xyz_j)
    apf.getParam(mesh.m_ptr, e, xi_j)
    me = apf.toModel(mesh.m_ptr, e)
    gmi.geval(g, me, xi_j, xyz2_j)

    for k=1:mesh.dim
      @test abs(xyz2_j[k] - xyz_j[k]) < 1e-5  # closestPoint is not that accurate
    end
  end

  return nothing
end

function test_update_coords(mesh, sbp, opts)

  @assert mesh.geoNums.can_eval

  # test updateCoords works
  xivec = getXiCoords(mesh)
  xvec = getXCoords(mesh)
  println("initially, xvec = \n", xvec)

  xivec .+= 1e-3*rand(length(xivec))
  update_coordsXi(mesh, sbp, opts, xivec)
  update_coords(mesh, sbp, opts, xvec)
  xvec2 = getXCoords(mesh)

  println("xvec = \n", xvec)
  println("xvec2 = \n", xvec2)
  println("diff = \n", xvec2 - xvec)
  @test maximum(abs.(xvec2 - xvec)) < 1e-12

  return nothing
end

function test_coords_file(opts, sbp, sbpface, mesh_ctor=PumiMeshDG2)

  mesh = mesh_ctor(Float64, sbp, opts, sbpface)

  # perturb coordinates
  if mesh.geoNums.can_eval
    xivec = getXiCoords(mesh)
    pert = 0.001*rand(size(xivec))
    xivec .+= pert
    update_coordsXi(mesh, sbp, opts, xivec)
  else
    xvec = getXCoords(mesh)
    pert = 0.01*rand(size(xvec))
    xvec .+= pert
    update_coords(mesh, sbp, opts, xvec)
  end

  writeCoordsBinary(mesh, "coordsb.dat")

  # load new mesh and check coordinates
  mesh2 = mesh_ctor(Float64, sbp, opts, sbpface)
  readCoordsBinary(mesh2, sbp, opts, "coordsb.dat")

  println("coords maxdiff = ", maximum(abs.(mesh.coords - mesh2.coords)))
  @test maximum(abs.(mesh.coords - mesh2.coords)) < 1e-15

  if mesh.geoNums.can_eval
    xivec2 = getXiCoords(mesh2)
    @test maximum(abs.(xivec2 - xivec)) < 1e-15
  end

  # test loading when mesh does not support xi coordinates
  opts2 = copy(opts)
  opts2["dmg_name"] = ".null"
  mesh3 = mesh_ctor(Float64, sbp, opts, sbpface)
  readCoordsBinary(mesh3, sbp, opts, "coordsb.dat")
  @test maximum(abs.(mesh.coords - mesh3.coords)) < 1e-15

  if mesh.geoNums.can_eval
    # save from mesh that does not support xi coordinates and load on mesh
    # that does
    writeCoordsBinary(mesh3, "coordsb2.dat")

    mesh4 = mesh_ctor(Float64, sbp, opts, sbpface)
    readCoordsBinary(mesh4, sbp, opts, "coordsb2.dat")

    @test maximum(abs.(mesh.coords - mesh4.coords)) < 1e-15

    xivec2 = getXiCoords(mesh4)
    @test maximum(abs.(xivec2 - xivec)) < 1e-7
  end

  #TODO: test Complex128
  return nothing
end

function test_refcounting()

  order = 1
  smb_name = "vortex.smb"
  dmg_name = "vortex.dmg"

  opts = Dict{Any, Any}(
    "dmg_name" => dmg_name,
    "smb_name" => smb_name,
    "order" => order,
    "coloring_distance" => 2,
    "numBC" => 1,
    "BC1" =>  [4, 7, 10, 13],
    "run_type" => 4,
    )
#    interp_op = [0.5 0 0; 0 0.5 0; 0 0 0.5]

  sbp = getTriSBPOmega(degree=order)
  vtx = sbp.vtx
  sbpface = TriFace{Float64}(order, sbp.cub, vtx)

  mesh = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=4)

  #fshape_ptr = apf.getConstantShapePtr(0)
  fshape_ptr = mesh.coordshape_ptr
  f_ptr = apf.createPackedField(mesh, "testfield1", 1, fshape_ptr)

  @test PdePumiInterface.countFieldRef(f_ptr) == 1

  f_ptr2 = apf.createPackedField(mesh, "testfield1", 1, fshape_ptr)
  @test f_ptr2 != f_ptr
  @test PdePumiInterface.countFieldRef(f_ptr2) == 1

  # load the same mesh from disk again, check the fields are not shared
  mesh2 = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=4)

  f_ptr3 = apf.createPackedField(mesh2, "testfield1", 1, fshape_ptr)

  @test mesh.m_ptr != mesh2.m_ptr
  for f in mesh.fields.orig
    @test !(f in mesh2.fields.orig)
  end
  for f in mesh.fields.user
    @test !(f in mesh2.fields.user)
  end

  m1_ptr = mesh.m_ptr
  finalize(mesh)
  @test apf.countMeshRefs(m1_ptr) == 0
  @test PdePumiInterface.countFieldRef(f_ptr) == 0

  # make sure finalizer is freeing fields
  finalize(mesh2)
  @test PdePumiInterface.countFieldRef(f_ptr3) == 0

  # create 2 Julia meshes sharing the same apf::Mesh
  mesh1 = PumiMeshDG2(Float64, sbp, opts, sbpface, dofpernode=4)

  #fshape_ptr = apf.getConstantShapePtr(0)
  fshape_ptr = mesh.coordshape_ptr
  f_orig_mesh1 = apf.createPackedField(mesh1.m_ptr, "origfield_mesh1",
                                       1, fshape_ptr)
  PdePumiInterface.attachOrigField(mesh1, f_orig_mesh1)

  n_orig_mesh1 = apf.createNumberingJ(mesh1.m_ptr, "orignum_mesh1", fshape_ptr,
                                      1)
  PdePumiInterface.attachOrigNumbering(mesh1, n_orig_mesh1)

  f_user_mesh1 = apf.createPackedField(mesh1.m_ptr, "userfield_mesh1",
                                       1, fshape_ptr)
  PdePumiInterface.attachUserField(mesh1, f_user_mesh1)

  n_user_mesh1 = apf.createNumberingJ(mesh1.m_ptr, "usernum_mesh1", fshape_ptr,
                                      1)
  PdePumiInterface.attachUserNumbering(mesh1, n_user_mesh1)

  mesh2 = PumiMeshDG2(mesh1, sbp, opts)

  n_user_mesh2 = apf.createNumberingJ(mesh2, "usernum_mesh2", 1, fshape_ptr)
  f_user_mesh2 = apf.createPackedField(mesh2, "userfield_mesh2", 1, fshape_ptr)

  # test orig fields are attached to both, but user fields are not
  @test f_orig_mesh1 in mesh2.fields.orig
  @test !(f_user_mesh1 in mesh2.fields.user)
  @test f_user_mesh1 in mesh1.fields.user

  @test n_orig_mesh1 in mesh2.numberings.orig
  @test !(n_user_mesh1 in mesh2.numberings.user)
  @test n_user_mesh1 in mesh1.numberings.user

  @test PdePumiInterface.countFieldRef(f_orig_mesh1) == 2
  @test PdePumiInterface.countFieldRef(f_user_mesh1) == 1

  @test PdePumiInterface.countNumberingRef(n_orig_mesh1) == 2
  @test PdePumiInterface.countNumberingRef(n_user_mesh1) == 1

  # write values to field
  vals = Float64[1.0]
  for (e, dim) in apf.FieldEntityIt(mesh1.m_ptr, fshape_ptr)
    apf.numberJ(n_orig_mesh1, e, 0, 0, 1)
    apf.numberJ(n_user_mesh1, e, 0, 0, 1)
    apf.numberJ(n_user_mesh2, e, 0, 0, 1)

    apf.setComponents(f_orig_mesh1, e, 0, vals)
    apf.setComponents(f_user_mesh1, e, 0, vals)
    apf.setComponents(f_user_mesh2, e, 0, vals)
  end

  writeVisFiles(mesh1, "mesh1_vis")
  writeVisFiles(mesh2, "mesh2_vis")

  mesh1_has_f_orig1 = false
  mesh1_has_n_orig1 = false
  mesh1_has_f_user1 = false
  mesh1_has_n_user1 = false

  mesh2_has_f_orig1 = false
  mesh2_has_n_orig1 = false
  mesh2_has_f_user1 = false
  mesh2_has_n_user1 = false

  mesh1_has_f_user2 = false
  mesh1_has_n_user2 = false
  mesh2_has_f_user2 = false
  mesh2_has_n_user2 = false

  for line in eachline("mesh1_vis/mesh1_vis.pvtu")
    if contains(line, "origfield_mesh1")
      mesh1_has_f_orig1 = true
    end
    if contains(line, "orignum_mesh1")
      mesh1_has_n_orig1 = true
    end
    if contains(line, "userfield_mesh1")
      mesh1_has_f_user1 = true
    end
    if contains(line, "usernum_mesh1")
      mesh1_has_n_user1 = true
    end
    if contains(line, "userfield_mesh2")
      mesh1_has_f_user2 = true
    end
    if contains(line, "usernum_mesh2")
      mesh1_has_n_user2 = true
    end
  end

  @test mesh1_has_f_orig1
  @test mesh1_has_n_orig1
  @test mesh1_has_f_user1
  @test mesh1_has_n_user1
  @test !mesh1_has_f_user2
  @test !mesh1_has_n_user2

  for line in eachline("mesh2_vis/mesh2_vis.pvtu")
    if contains(line, "origfield_mesh1")
      mesh2_has_f_orig1 = true
    end
    if contains(line, "orignum_mesh1")
      mesh2_has_n_orig1 = true
    end
    if contains(line, "userfield_mesh1")
      mesh2_has_f_user1 = true
    end
    if contains(line, "usernum_mesh1")
      mesh2_has_n_user1 = true
    end
    if contains(line, "userfield_mesh2")
      mesh2_has_f_user2 = true
    end
    if contains(line, "usernum_mesh2")
      mesh2_has_n_user2 = true
    end
  end

  @test mesh2_has_f_orig1
  @test mesh2_has_n_orig1
  @test !mesh2_has_f_user1
  @test !mesh2_has_n_user1
  @test mesh2_has_f_user2
  @test mesh2_has_n_user2

  finalize(mesh1)
  @test PdePumiInterface.countFieldRef(f_orig_mesh1) == 1
  @test PdePumiInterface.countFieldRef(f_user_mesh1) == 0

  @test PdePumiInterface.countNumberingRef(n_orig_mesh1) == 1
  @test PdePumiInterface.countNumberingRef(n_user_mesh1) == 0

  finalize(mesh2)

  return nothing
end

function test_interioredge(mesh)
# this only works with the mesh corner20_split_mesh0_.smb

  nfaces = 12 +12 + 9 + 8 + 4
  @test mesh.numBoundaryFaces == nfaces
  @test length(mesh.bndryfaces) == nfaces

  for face in mesh.bndryfaces
    @test face.element > 0
    @test face.face <= mesh.numEl
    @test face.face > 0
    @test face.face <= 3
  end

  # test no duplicates
  bndryfaces = sort(mesh.bndryfaces)
  for i=2:mesh.numBoundaryFaces
    @test bndryfaces[i] != bndryfaces[i - 1]
  end

  return nothing
end

function test_split(mesh::PumiMesh)
  # test splitting a mesh

  nshape = apf.getConstantShapePtr(mesh.dim)
  partnums = apf.createNumberingJ(mesh, "split_partnums", 1, nshape)

  apf.getDefaultSplit(mesh.m_ptr, 4, partnums)

  # count the number of element assigned to each part, make sure it satisfies
  # some tolerance
  elcounts = zeros(Int, 4)
  for el in mesh.elements
    partnum = apf.getNumberJ(partnums, el, 0, 0) + 1
    elcounts[partnum] += 1
  end

  println("elcounts = ", elcounts)
  expected_count = mesh.numEl/4
  @test sum(elcounts) == mesh.numEl
  for i=1:4
    @test elcounts[i] >= expected_count/1.1
    @test elcounts[i] <= 1.1*expected_count
  end

  test_split_submesh(mesh)
  return nothing
end

function test_split_submesh(mesh::PumiMesh)
  # test creating a submesh and splitting it

  println("testing splitting submesh")
  # create submesh with first half of elements
  nel = div(mesh.numEl, 2)
  elnums = zeros(Cint, nel)
  for i=1:nel
    elnums[i] = i
  end
  sdata = apf.createSubMesh(mesh.m_ptr, mesh.entity_Nptrs, elnums)
  m2_ptr = apf.getNewMesh(sdata)

  # create the apf::Numbering on the submesh
  nshape = apf.getConstantShapePtr(mesh.dim)
  partnums = apf.createNumberingJ(m2_ptr, "split_partnums2", nshape, 1)

  apf.getDefaultSplit(m2_ptr, 4, partnums)

  # count the number of element assigned to each part, make sure it satisfies
  # some tolerance
  elcounts = zeros(Int, 4)
  for el in apf.MeshIterator(m2_ptr, mesh.dim)
    partnum = apf.getNumberJ(partnums, el, 0, 0) + 1
    elcounts[partnum] += 1
  end

  expected_count = nel/4
  @test sum(elcounts) == nel
  for i=1:4
    @test elcounts[i] >= expected_count/1.1
    @test elcounts[i] <= 1.1*expected_count
  end

  apf.free(sdata)

  return nothing
end

function test_submesh_transfer(mesh::PumiMesh)
# test transferring numberings to submesh

  nel = div(mesh.numEl, 2)
  elnums = zeros(Cint, nel)
  for i=1:nel
    elnums[i] = i
  end
  sdata = apf.createSubMesh(mesh.m_ptr, mesh.entity_Nptrs, elnums)
  m2_ptr = apf.getNewMesh(sdata)

  n_old_ptr = mesh.vert_Nptr
  n_new_ptr = apf.transferNumbering(sdata, n_old_ptr)

  # verify that vertices with same coords have same numbers
  coords_old = zeros(Float64, 3)
  coords_new = zeros(Float64, 3)
  for e_old in apf.MeshIterator(mesh.m_ptr, 0)
    apf.getPoint(mesh.m_ptr, e_old, 0, coords_old)
    for e_new in apf.MeshIterator(m2_ptr, 0)
      apf.getPoint(m2_ptr, e_new, 0, coords_new)

      if norm(coords_old - coords_new) < 1e-10
        num_old = apf.getNumberJ(n_old_ptr, e_old, 0, 0)
        num_new = apf.getNumberJ(n_new_ptr, e_new, 0, 0)
        @test num_old == num_new
      end
    end
  end

  apf.free(sdata)

  return nothing
end


text:
function nonsym_xxz_ham(;spin = 1,delta = 1,zfield = 0.0)
    (sx,sy,sz,id) = nonsym_spintensors(spin)

    mpo=Array{typeof(sx),3}(undef,1,5,5)
    mpo[1,:,:]=[id sx sy delta*sz zfield*sz;
                0*id 0*id 0*id 0*id sx;
                0*id 0*id 0*id 0*id sy;
                0*id 0*id 0*id 0*id sz;
                0*id 0*id 0*id 0*id id]

    return MPOHamiltonian(mpo)
end

function su2_xxx_ham(;spin = 1//2)
    #only checked for spin = 1 and spin = 2...
    ph = ℂ[SU₂](spin=>1)

    Sl1 = TensorMap(ones, Defaults.eltype, ℂ[SU₂](0=>1)*ph , ℂ[SU₂](1=>1)*ph)*sqrt(spin^2+spin)
    Sr1 = TensorMap(ones, Defaults.eltype, ℂ[SU₂](1=>1)*ph , ℂ[SU₂](0=>1)*ph)*sqrt(spin^2+spin)

    return MPOHamiltonian([Sl1,Sr1]);
end

function u1_xxz_ham(;spin = 1,delta = 1,zfield = 0.0)
    (sxd,syd,szd,idd) = spinmatrices(spin);
    @tensor ham[-1 -2;-3 -4]:=sxd[-1,-3]*sxd[-2,-4]+syd[-1,-3]*syd[-2,-4]+(delta*szd)[-1,-3]*szd[-2,-4]+zfield*0.5*szd[-1,-3]*idd[-2,-4]+zfield*0.5*idd[-1,-3]*szd[-2,-4]

    indu1map = [U₁(v) for v in -spin:1:spin];
    pspace = U1Space((v=>1 for v in indu1map));

    symham = TensorMap(zeros,eltype(ham),pspace*pspace,pspace*pspace)

    for (i,j,k,l) in Iterators.product(1:size(ham,1),1:size(ham,1),1:size(ham,1),1:size(ham,1))
        if ham[i,j,k,l]!=0
            copyto!(symham[(indu1map[i],indu1map[j],indu1map[k],indu1map[l])],ham[i,j,k,l])
        end
    end

    return MPOHamiltonian(decompose_localmpo(add_util_leg(symham)))
end

ext: jl


text:
using TableWidgets, Observables, WebIO, Widgets, InteractBase
using DataFrames
using Test

@testset "selector" begin
    v = [1, 1, 2, 2, 1, 3, 4]
    sel = categoricalselector(v)
    @test observe(sel)[] == v
    observe(sel, :checkboxes)[] = [1, 3]
    sleep(0.1)
    @test observe(sel)[] == [1, 1, 1, 3]

    sel = rangeselector(v)
    @test observe(sel)[] == v
    observe(sel, :extrema)[] = [1.8, 3.2]
    observe(sel, :changes)[] = 4
    sleep(0.1)
    @test observe(sel)[] == [2, 2, 3]

    sel = selector(v, map)
    @test observe(sel)[] == fill(true, length(v))
    observe(sel, :function)[] = t -> t != 2
    observe(sel, :textbox, :changes)[] = 4
    sleep(0.1)
    @test observe(sel)[] == map(t -> t != 2, v)
end

@testset "undo" begin
    obs = Observable(1)
    u = TableWidgets.Undo(obs)
    obs[] = 2
    sleep(0.1)
    @test obs[] == 2
    u()
    sleep(0.1)
    @test obs[] == 1
    @test observe(u)[] == 1
end

@testset "table" begin
    df = DataFrame(x = 1:4, y = ["a", "b", "c", "d"])
    n = Observable(3)
    wdg = TableWidgets.head(df, n)
    l = WebIO.children(WebIO.children(wdg[:head][].dom)[2]) |> length
    @test l == 3
    n[] = 10
    sleep(0.1)
    l = WebIO.children(WebIO.children(wdg[:head][].dom)[2]) |> length
    @test l == 4
end


text:
using Test, MPI
import GaussQuadrature
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.Mesh.Geometry
using CLIMA.Mesh.Interpolation
using StaticArrays

#------------------------------------------------
using CLIMA.DGmethods
using CLIMA.DGmethods.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.LowStorageRungeKuttaMethod
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using CLIMA.Atmos
using CLIMA.VariableTemplates
using CLIMA.MoistThermodynamics
using CLIMA.PlanetParameters
using CLIMA.TicToc
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates
using CLIMA.VTK

using CLIMA.Atmos: vars_state, vars_aux

using Random
using Statistics
const seed = MersenneTwister(0)

const ArrayType = CLIMA.array_type()

#------------------------------------------------
if !@isdefined integration_testing
    const integration_testing =
    parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_INTEGRATION_TESTING","false")))
end
#------------------------------------------------
function run_brick_interpolation_test()
  MPI.Initialized() || MPI.Init()

#@testset "LocalGeometry" begin
    FT = Float64
    ArrayType = Array
    mpicomm = MPI.COMM_WORLD

    xmin, ymin, zmin = 0, 0, 0                   # defining domain extent
    xmax, ymax, zmax = 2000, 400, 2000
    xres = [FT(200), FT(200), FT(200)] # resolution of interpolation grid

    xgrd = range(xmin, xmax, step=xres[1])
    ygrd = range(ymin, ymax, step=xres[2])
    zgrd = range(zmin, zmax, step=xres[3])

 #   Ne        = (20,2,20)
    Ne        = (4,2,4)

    polynomialorder = 8 #8 #4
    #-------------------------
    _x, _y, _z = CLIMA.Mesh.Grids.vgeoid.x1id, CLIMA.Mesh.Grids.vgeoid.x2id, CLIMA.Mesh.Grids.vgeoid.x3id
    _ρ, _ρu, _ρv, _ρw = 1, 2, 3, 4
    #-------------------------

    brickrange = (range(FT(xmin); length=Ne[1]+1, stop=xmax),
                  range(FT(ymin); length=Ne[2]+1, stop=ymax),
                  range(FT(zmin); length=Ne[3]+1, stop=zmax))
    topl = StackedBrickTopology(MPI.COMM_SELF, brickrange, periodicity = (true, true, false))

    grid = DiscontinuousSpectralElementGrid(topl,
                                            FloatType = FT,
                                            DeviceArray = ArrayType,
                                            polynomialorder = polynomialorder)

    model = AtmosModel(FlatOrientation(),
                     NoReferenceState(),
					 ConstantViscosityWithDivergence(FT(0)),
                     EquilMoist(),
                     NoRadiation(),
                     NoSubsidence{FT}(),
                     (Gravity()),
					 NoFluxBC(),
                     Initialize_Brick_Interpolation_Test!)

    dg = DGModel(model,
               grid,
               Rusanov(),
               CentralNumericalFluxDiffusive(),
               CentralGradPenalty())

    Q = init_ode_state(dg, FT(0))
    #------------------------------
    x1 = @view grid.vgeo[:,_x,:]
    x2 = @view grid.vgeo[:,_y,:]
    x3 = @view grid.vgeo[:,_z,:]

    st_idx = _ρ # state vector
    elno = 10

    var = @view Q.data[:,st_idx,:]

    #fcn(x,y,z) = x .* y .* z # sample function
    fcn(x,y,z) = sin.(x) .* cos.(y) .* cos.(z) # sample function

    var .= fcn( x1 ./ xmax, x2 ./ ymax, x3 ./ zmax )
    #----calling interpolation function on state variable # st_idx--------------------------
    #intrp_brck = InterpolationBrick(grid, xres, FT)
    intrp_brck = InterpolationBrick(grid, xres)
    interpolate_brick!(intrp_brck, Q.data, st_idx)

    #------testing
    Nel = length( grid.topology.realelems )

    error = zeros(FT, Nel)

    for elno in 1:Nel
      fex = similar(intrp_brck.V[elno])
      fex = fcn( intrp_brck.x[elno][1,:] ./ xmax , intrp_brck.x[elno][2,:] ./ ymax , intrp_brck.x[elno][3,:] ./ zmax )
      error[elno] = maximum(abs.(intrp_brck.V[elno][:]-fex[:]))
    end

    l_infinity_local = maximum(error)
    l_infinity_domain = MPI.Allreduce(l_infinity_local, MPI.MAX, mpicomm)

    return l_infinity_domain < 1.0e-14
    #----------------
end #function run_brick_interpolation_test

#-----taken from Test example
function Initialize_Brick_Interpolation_Test!(state::Vars, aux::Vars, (x,y,z), t)
    FT         = eltype(state)

    # Dummy variables for initial condition function
    state.ρ     = FT(0)
    state.ρu    = SVector{3,FT}(0,0,0)
    state.ρe    = FT(0)
    state.moisture.ρq_tot = FT(0)
end
#------------------------------------------------
Base.@kwdef struct TestSphereSetup{FT}
  p_ground::FT = MSLP
  T_initial::FT = 255
  domain_height::FT = 30e3
end
#----------------------------------------------------------------------------
# Cubed sphere, lat/long interpolation test
#----------------------------------------------------------------------------
function run_cubed_sphere_interpolation_test()
    CLIMA.init()

    FT = Float64
    mpicomm = MPI.COMM_WORLD
    root = 0

    ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
    loglevel = Dict("DEBUG" => Logging.Debug,
                    "WARN"  => Logging.Warn,
                    "ERROR" => Logging.Error,
                    "INFO"  => Logging.Info)[ll]

    logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
    global_logger(ConsoleLogger(logger_stream, loglevel))

    domain_height = FT(30e3)

    polynomialorder = 12#1#4 #5
    numelem_horz = 3#4 #6
    numelem_vert = 4#1 #1 #1#6 #8

    #-------------------------
    _x, _y, _z = CLIMA.Mesh.Grids.vgeoid.x1id, CLIMA.Mesh.Grids.vgeoid.x2id, CLIMA.Mesh.Grids.vgeoid.x3id
    _ρ, _ρu, _ρv, _ρw = 1, 2, 3, 4
    #-------------------------
    vert_range = grid1d(FT(planet_radius), FT(planet_radius + domain_height), nelem = numelem_vert)

 #   vert_range = grid1d(FT(1.0), FT(2.0), nelem = numelem_vert)

    lat_res  = 5 * π / 180.0 # 5 degree resolution
    long_res = 5 * π / 180.0 # 5 degree resolution
    r_res    = (vert_range[end] - vert_range[1])/FT(numelem_vert) #1000.00    # 1000 m vertical resolution

    #----------------------------------------------------------
    setup = TestSphereSetup{FT}()

    topology = StackedCubedSphereTopology(mpicomm, numelem_horz, vert_range)

    grid = DiscontinuousSpectralElementGrid(topology,
                                            FloatType = FT,
                                            DeviceArray = ArrayType,
                                            polynomialorder = polynomialorder,
                                            meshwarp = CLIMA.Mesh.Topologies.cubedshellwarp)

    model = AtmosModel(SphericalOrientation(),
                       NoReferenceState(),
                       ConstantViscosityWithDivergence(FT(0)),
                       DryModel(),
                       NoRadiation(),
                       NoSubsidence{FT}(),
                       nothing,
                       NoFluxBC(),
                       setup)

    dg = DGModel(model, grid, Rusanov(),
                 CentralNumericalFluxDiffusive(), CentralGradPenalty())

    Q = init_ode_state(dg, FT(0))
    #------------------------------
    x1 = @view grid.vgeo[:,_x,:]
    x2 = @view grid.vgeo[:,_y,:]
    x3 = @view grid.vgeo[:,_z,:]

    xmax = maximum( abs.(x1) )
    ymax = maximum( abs.(x2) )
    zmax = maximum( abs.(x3) )

    st_idx = _ρ # state vector

    var = @view Q.data[:,st_idx,:]

#    fcn(x,y,z) = x .* y .* z # sample function
    fcn(x,y,z) = sin.(x) .* cos.(y) .* cos.(z) # sample function

    var .= fcn( x1 ./ xmax, x2 ./ ymax, x3 ./ zmax )
  #------------------------------
    @time intrp_cs = InterpolationCubedSphere(grid, collect(vert_range), numelem_horz, lat_res, long_res, r_res)
    interpolate_cubed_sphere!(intrp_cs, Q.data, st_idx)
    #----------------------------------------------------------

    Nel = length( grid.topology.realelems )

    error = zeros(FT, Nel)

    for elno in 1:Nel
        if ( length(intrp_cs.V[elno]) > 0 )
            fex = similar(intrp_cs.V[elno])
            x1g = similar(intrp_cs.V[elno])
            x2g = similar(intrp_cs.V[elno])
            x3g = similar(intrp_cs.V[elno])

            x1_grd = intrp_cs.radc[elno] .* sin.(intrp_cs.latc[elno]) .* cos.(intrp_cs.longc[elno]) # inclination -> latitude; azimuthal -> longitude.
            x2_grd = intrp_cs.radc[elno] .* sin.(intrp_cs.latc[elno]) .* sin.(intrp_cs.longc[elno]) # inclination -> latitude; azimuthal -> longitude.
            x3_grd = intrp_cs.radc[elno] .* cos.(intrp_cs.latc[elno])

            fex = fcn( x1_grd ./ xmax , x2_grd ./ ymax , x3_grd ./ zmax )
            error[elno] = maximum(abs.(intrp_cs.V[elno][:]-fex[:]))
        end
    end
    #----------------------------------------------------------
    l_infinity_local = maximum(error)
    l_infinity_domain = MPI.Allreduce(l_infinity_local, MPI.MAX, mpicomm)

#----------------------------------------------------------------------------
    return l_infinity_domain < 1.0e-12
end
#----------------------------------------------------------------------------
function (setup::TestSphereSetup)(state, aux, coords, t)
  # callable to set initial conditions
  FT = eltype(state)

  r = norm(coords, 2)
  h = r - FT(planet_radius)

  scale_height = R_d * setup.T_initial / grav
  p = setup.p_ground * exp(-h / scale_height)

  state.ρ = air_density(setup.T_initial, p)
  state.ρu = SVector{3, FT}(0, 0, 0)
  state.ρe = state.ρ * (internal_energy(setup.T_initial) + aux.orientation.Φ)
  nothing
end
#----------------------------------------------------------------------------

@testset "Interpolation tests" begin
    @test run_brick_interpolation_test()
    @test run_cubed_sphere_interpolation_test()
end
#------------------------------------------------

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
#==
Run example with default log level:

    >julia simple.jl

Run example with elevated log level:

    >LOG4JL_DEFAULT_STATUS_LEVEL=ALL julia simple.jl

==#
module X
    println("Current:", current_module())
    # println("Current:", module_parent(current_module()))

    using Log4jl
    @Log4jl logger = @logger

    module Y
        println("Current:", current_module())

        using Log4jl
        @Log4jl logger = @logger(URI="log4jl.json")
        Log4jl.level!(logger, Log4jl.Level.WARN)

        module Z
            println("Current:", current_module())

            using Log4jl
            @Log4jl logger = @Log4jl.logger begin
                # Create a configuration for the logger
                cfg = Log4jl.DefaultConfiguration()
                # Add an appender to the configuration
                Log4jl.appender!(cfg, OUTPUT = Log4jl.Appenders.Console(
                    layout = Log4jl.Layouts.PatternLayout("%d{%Y-%m-%d %H:%M:%S} [%t] %-5p %l %c{3} - %m%n")
                ))
                # Add logger configuration with the name of the logger (i.e. "X.Y.Z")
                #Log4jl.logger!(cfg, "X.Y.Z", Log4jl.LoggerConfig("LWC", Log4jl.Level.WARN))
                Log4jl.logger!(cfg, "X.Y.Z", Log4jl.Level.WARN)
                # Reference the appender in the logger configuration
                Log4jl.reference!(cfg, "X.Y.Z" => "OUTPUT")
                return cfg
            end

            println(logger)
            println(logger.config)

            function test_log3()
                fatal(logger, "L3: fatal {}", 5)
                @fatal "L3: FATAL {}" 5
                @debug "L3: DEBUG {}" 1
                @warn  "L3: WARN {}" 3
                @info  "L3: INFO {}" 2
                @debug "L3: DEBUG"
                trace(logger, "L3: TRACE")
            end

        end

        function test_log2()
            Y.Z.test_log3()

            @trace "L2: TRACE"
            @debug "L2: DEBUG"
            @info  "L2: INFO {}" 2
            @warn  "L2: WARN {}" 3
            @error "L2: ERROR {}" 4
            @fatal "L2: FATAL {}" 5
        end

    end

    function test_log1()
        Y.test_log2()

        trace(logger, "aaaaa")
        debug(logger, "bbbbb {}", 1)
        info(logger,  "ccccc {}", 2)
        @warn "ddddd {}" 3
        sleep(1.)
        fatal(logger, "fffff {}", 5)
        error(logger, "eeeee {}", 4)
        Log4jl.verbose(logger, "vvvvv {}", 6)
    end

end

println("Current:", current_module())

println("Registered contexts:")
for (ctxname, ctx) in Log4jl.contexts(Log4jl.LOG4JL_CONTEXT_SELECTOR)
    println("\t$ctxname => $ctx")
    for (lgrname, lgr) in ctx.loggers
        # lgr2 = Log4jl.logger(ctx, lgrname)
        println("\t\t$lgrname => $lgr")
        println("\t\t\t$(lgr.config) => $(Log4jl.references(lgr.config))")
    end
end
println()

X.test_log1()
println("END\n")

xdump(X.logger)


text:
tests = [
    "planetoid",
    "cora",
    "ppi",
    "reddit",
    "qm7b",
    "entities"
]

@testset "datasets" begin
    @test_throw ArgumentError traindata(Reddit())
    @test_throw ArgumentError validdata(Reddit())
    @test_throw ArgumentError testdata(Reddit())
    @test_throw AssertionError traindata(Planetoid(), :abc)

    for t in tests
        include("$(t).jl")
    end
end


text:
#=

# 204: 2D Convection in Hagen-Poiseuille flow
([source code](SOURCE_URL))

Solve the equation

```math
\partial_t u -\nabla ( D \nabla u - v u) = 0
```
in $\Omega=(0,L)\times (0,H)$ with dirichlet boundary conditions
at $x=0$ and outflow boundary condition at $x=L$.
=#

module Example204_HagenPoiseuille
using Printf
using VoronoiFVM
using ExtendableGrids
using GridVisualize

## Bernoulli function used in the exponential fitting discretization
function bernoulli(x)
    if abs(x)<nextfloat(eps(typeof(x)))
        return 1
    end
    return x/(exp(x)-1)
end

function exponential_flux!(f,u0,edge,data)
    u=unknowns(edge,u0)
    vh=project(edge,data.v)
    Bplus= data.D*bernoulli(vh/data.D)
    Bminus=data.D*bernoulli(-vh/data.D)
    f[1]=Bminus*u[1,1]-Bplus*u[1,2]
end

function outflow!(f,u,node,data)
    if node.region==2
        f[1]=data.v[1]*u[1]
    end
end

function main(;nref=0,Plotter=nothing,D=0.01,v=1.0,tend=100,cin=1.0)
    H=1.0
    L=5.0
    grid=simplexgrid(range(0,L,length=20*2^nref),
                     range(0,H,length=5*2^nref))

    function fhp(x,y)
        yh=y/H
        return v*4*yh*(1.0-yh),0
    end

    function flux!(f,u,edge)
        vd=evelo[edge.index]/D
        bp=fbernoulli(vd)
        bm=fbernoulli(-vd)
        f[1]=D*(bp*u[1] - bm*u[2])
    end

    function outflow!(f,u,node)
        if node.region==2
            f[1]=bfvelo[node.ibnode,node.ibface]*u[1]
        end
    end

    ispec=1
    physics=VoronoiFVM.Physics(num_species=1,flux=flux!,breaction=outflow!)
    sys=VoronoiFVM.DenseSystem(grid,physics)
    enable_species!(sys,ispec,[1])
    evelo=edgevelocities(grid,fhp)
    bfvelo=bfacevelocities(grid,fhp)

    boundary_dirichlet!(sys,ispec,4,cin)

    ## Create a solution array
    inival=unknowns(sys,inival=0)

    ## Transient solution of the problem
    control=VoronoiFVM.NewtonControl()
    control.Δt=0.01*2.0^(-nref)
    control.Δt_min=0.01*2.0^(-nref)
    control.Δt_max=0.1*tend
    tsol=solve(inival,sys,[0,tend],control=control)

    vis=GridVisualizer(Plotter=Plotter)
    for i=1:length(tsol.t)
        scalarplot!(vis[1,1],grid,tsol[1,:,i],flimits=(0,cin+1.0e-5))
    end
    tsol
end

function test()
    tsol=main()
    all(tsol[end].≈1)
end

end


text:
module AerostructuralDynamics

using ArnoldiMethod
using ForwardDiff
using GXBeam
using LinearAlgebra
using LinearMaps
using FillArrays
using StaticArrays
using DiffEqBase
using RecipesBase
using Printf

# Models
export Steady
export QuasiSteady
export Wagner
export Peters
export LiftingLine
export TypicalSection
export RigidBody
export GEBT
export SimpleFlap
export LiftingLineFlaps
export Trim

# Interface
export couple_models
export number_of_states
export number_of_inputs
export number_of_parameters
export state_indices
export input_indices
export parameter_indices
export get_states, set_states!
export get_inputs, set_inputs!
export get_parameters, set_parameters!
export get_additional_parameters, set_additional_parameters!
export separate_states
export separate_inputs
export separate_parameters
export get_residual, get_residual!
export get_rate_jacobian, get_rate_jacobian!
export get_state_jacobian, get_state_jacobian!
export get_input_jacobian, get_input_jacobian!
export get_parameter_jacobian, get_parameter_jacobian!
export get_time_gradient, get_time_gradient!
export get_coupling_inputs, get_coupling_inputs!
export get_eigen
export get_ode
export correlate_eigenmodes

# needed for traits
import Base.isempty
import Base.iszero

# Interface Functions
include("traits.jl")
include("interface.jl")
include("mode-tracking.jl")

# Two-Dimensional Models

# Aerodynamic Models
include(joinpath("aerodynamics", "quasisteady.jl"))
include(joinpath("aerodynamics", "wagner.jl"))
include(joinpath("aerodynamics", "peters.jl"))

# Structural Models
include(joinpath("structures", "section.jl"))
include(joinpath("structures", "liftingline-section.jl"))

# Control Surface Models
include(joinpath("control-surfaces", "simpleflap.jl"))

# Controller Models
include(joinpath("controllers", "liftingline-section-control.jl"))

# Coupled Models
include(joinpath("couplings", "two-dimensional", "quasisteady-section.jl"))
include(joinpath("couplings", "two-dimensional", "wagner-section.jl"))
include(joinpath("couplings", "two-dimensional", "peters-section.jl"))

include(joinpath("couplings", "two-dimensional", "quasisteady-section-simpleflap.jl"))
include(joinpath("couplings", "two-dimensional", "wagner-section-simpleflap.jl"))
include(joinpath("couplings", "two-dimensional", "peters-section-simpleflap.jl"))

include(joinpath("couplings", "two-dimensional", "quasisteady-liftingline.jl"))
include(joinpath("couplings", "two-dimensional", "wagner-liftingline.jl"))
include(joinpath("couplings", "two-dimensional", "peters-liftingline.jl"))

include(joinpath("couplings", "two-dimensional", "quasisteady-liftingline-simpleflap.jl"))
include(joinpath("couplings", "two-dimensional", "wagner-liftingline-simpleflap.jl"))
include(joinpath("couplings", "two-dimensional", "peters-liftingline-simpleflap.jl"))

# Three-Dimensional Models

# Aerodynamic Models
include(joinpath("aerodynamics", "liftingline.jl"))

# Structural Models
include(joinpath("structures", "gxbeam.jl"))

# Dynamics Models
include(joinpath("dynamics", "rigidbody.jl"))

# Control Surface Models
include(joinpath("control-surfaces", "liftinglineflaps.jl"))

# Controller Models
include(joinpath("controllers", "trim.jl"))

# Coupled Models
include(joinpath("couplings", "three-dimensional", "liftingline-rigidbody.jl"))
include(joinpath("couplings", "three-dimensional", "liftingline-gxbeam.jl"))
include(joinpath("couplings", "three-dimensional", "liftingline-gxbeam-rigidbody.jl"))

include(joinpath("couplings", "three-dimensional", "liftingline-rigidbody-liftinglineflaps.jl"))
include(joinpath("couplings", "three-dimensional", "liftingline-gxbeam-liftinglineflaps.jl"))
include(joinpath("couplings", "three-dimensional", "liftingline-gxbeam-rigidbody-liftinglineflaps.jl"))

include(joinpath("couplings", "three-dimensional", "liftingline-rigidbody-liftinglineflaps-trim.jl"))
include(joinpath("couplings", "three-dimensional", "liftingline-gxbeam-liftinglineflaps-trim.jl"))
include(joinpath("couplings", "three-dimensional", "liftingline-gxbeam-rigidbody-liftinglineflaps-trim.jl"))

end


text:
using Turing: Tracker

function turing_inference(
    prob::DiffEqBase.DEProblem,
    alg,
    t,
    data,
    priors;
    likelihood_dist_priors = [InverseGamma(2, 3)],
    likelihood = (u,p,t,σ) -> MvNormal(u, σ[1]*ones(length(u))),
    num_samples=1000, sampler = Turing.NUTS(0.65),
    syms = [Turing.@varname(theta[i]) for i in 1:length(priors)],
    sample_u0 = false,
    progress = true,
    solve_progress = false,
    cb = θ -> false,
    kwargs...,
)
    N = length(priors)
    Turing.@model mf(x, ::Type{T} = Float64) where {T <: Real} = begin
        theta = Vector{T}(undef, length(priors))
        for i in 1:length(priors)
            theta[i] ~ NamedDist(priors[i], syms[i])
        end
        σ = Vector{T}(undef, length(likelihood_dist_priors))
        for i in 1:length(likelihood_dist_priors)
            σ[i] ~ likelihood_dist_priors[i]
        end
        nu = length(prob.u0)
        u0 = convert.(T, sample_u0 ? theta[1:nu] : prob.u0)
        p = convert.(T, sample_u0 ? theta[(nu + 1):end] : theta)
        _saveat = isnothing(t) ? Float64[] : t
        sol = concrete_solve(prob, alg, u0, p; saveat = _saveat, progress = solve_progress, kwargs...)
        failure = size(sol, 2) < length(_saveat)
        if failure
            @logpdf() = T(0) * sum(p) + T(-Inf)
            return
        end
        if ndims(sol) == 1
            x ~ likelihood(sol[:], theta, Inf, σ)
        else
            for i = 1:length(t)
                x[:, i] ~ likelihood(sol[:, i], theta, _saveat[i], σ)
            end
        end
        cb(Tracker.data.([u0; p; σ]))
        return
    end

    # Instantiate a Model object.
    model = mf(data)
    if num_samples > 0
        chn = sample(model, sampler, num_samples; progress = progress)
    else
        chn = Chains(zeros(0, 0, 0))
    end
    vi = Turing.VarInfo(model)
    backend = try
        Turing.Core.getADbackend(sampler)
    catch e
        Turing.ADBackend(:forwarddiff)()
    end
    function nlogp(θ)
        spl = Turing.SampleFromPrior()
        vi′ = Turing.VarInfo(vi, spl, θ)
        model(vi′, spl); cb(θ)
        -Turing.getlogp(vi′)
    end
    function ∇nlogp(θ)
        -Turing.gradient_logp(backend, θ, vi, model)[2]
    end
    return setinfo(chn, (model = model, nlogp = nlogp, ∇nlogp = ∇nlogp))
end


text:
struct Graph{N}
    directchildren::Vector{Vector{Int64}} # direct successor nodes
    loopchildren::Vector{Vector{Int64}} # successor nodes excluding direct successors
    ineqchildren::Vector{Vector{Int64}}
    successors::Vector{Vector{Int64}} # direct and loop children
    predecessors::Vector{Vector{Int64}} # direct parent and loop-opening predecessor(s?) (for numerics?)
    connections::Vector{Vector{Int64}} # direct connections

    dfslist::SVector{N,Int64} # depth-first-seach list (dfslist[end] = root)
    rdfslist::SVector{N,Int64} # reverse dfslist

    dict::UnitDict{Base.OneTo{Int64},Int64} # maps ids to the graph-interal numbering (not dfs order)
    rdict::UnitDict{Base.OneTo{Int64},Int64} # reverse mapping
    activedict::UnitDict{Base.OneTo{Int64},Bool}

    function Graph(origin::Origin,bodies::Vector{<:Body},
            eqconstraints::Vector{<:EqualityConstraint},ineqconstraints::Vector{<:InequalityConstraint}
        )

        oid = origin.id
        adjacency, dict = adjacencyMatrix(eqconstraints, bodies)
        dfsgraph, dfslist, loops = dfs(adjacency, dict, oid)
        pat = pattern(dfsgraph, dict, loops)
        fil, originals = fillins(dfsgraph, pat, dict, loops)

        adjacency = deleteat(adjacency, dict[oid])
        dfsgraph = deleteat(dfsgraph, dict[oid])
        pat = deleteat(pat, dict[oid])
        fil = deleteat(fil, dict[oid])
        originals = deleteat(originals, dict[oid])
        dfslist = StaticArrays.deleteat(dfslist, length(dfslist))

        for (id, ind) in dict
            ind > dict[oid] && (dict[id] = ind - 1)
        end
        pop!(dict, oid)
        rdict = Dict(ind => id for (id, ind) in dict)
        activedict = Dict(id => true for (id, ind) in dict)

        for eqc in eqconstraints
            eqc.parentid == oid && (eqc.parentid = nothing)
            activedict[eqc.id] = isactive(eqc)
        end

        for body in bodies
            activedict[body.id] = isactive(body)
        end

        for ineqc in ineqconstraints
            activedict[ineqc.id] = isactive(ineqc)
        end

        N = length(dict)

        adjacency = convert(Vector{SVector{N,Bool}}, adjacency)
        dfsgraph = convert(Vector{SVector{N,Bool}}, dfsgraph)
        pat = convert(Vector{SVector{N,Bool}}, pat)
        fil = convert(Vector{SVector{N,Bool}}, fil)
        originals = convert(Vector{SVector{N,Bool}}, originals)

        dirs = directchildren(dfslist, originals, dict)
        loos = loopchildren(dfslist, fil, dict)
        ineqs = ineqchildren(dfslist, bodies, ineqconstraints, dict)
        sucs = successors(dfslist, pat, dict)
        preds = predecessors(dfslist, pat, dict)
        cons = connections(dfslist, adjacency, dict)

        dict = UnitDict(dict)
        rdict = UnitDict(rdict)
        activedict = UnitDict(activedict)

        new{N}(dirs, loos, ineqs, sucs, preds, cons, dfslist, reverse(dfslist), dict, rdict, activedict)
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, graph::Graph{N}) where {N}
    summary(io, graph)
end

### Construction functions
function adjacencyMatrix(eqconstraints::Vector{<:EqualityConstraint}, bodies::Vector{<:Body})
    A = zeros(Bool, 0, 0)
    dict = Dict{Int64,Int64}()
    n = 0

    for constraint in eqconstraints
        childid = constraint.id
        A = [A zeros(Bool, n, 1); zeros(Bool, 1, n) zero(Bool)]
        dict[childid] = n += 1
        for bodyid in unique([constraint.parentid;constraint.childids])
            if !haskey(dict, bodyid)
                A = [A zeros(Bool, n, 1); zeros(Bool, 1, n) zero(Bool)]
                dict[bodyid] = n += 1
            end
            A[dict[childid],dict[bodyid]] = true
        end
    end
    for body in bodies # add unconnected bodies
        if !haskey(dict, body.id)
            A = [A zeros(Bool, n, 1); zeros(Bool, 1, n) zero(Bool)]
            dict[body.id] = n += 1
        end
    end

    A = A .| A'
    return convert(Matrix{Bool}, A), dict
end

function dfs(adjacency::Matrix, dict::Dict, originid::Integer)
    N = size(adjacency)[1]
    dfsgraph = zeros(Bool, N, N)
    dfslist = zeros(Int64, N)
    visited = zeros(Bool, N)
    loops = Vector{Int64}[]
    index = N

    dfslist[index] = originid
    visited[dict[originid]] = true
    dfs!(adjacency, dfsgraph, dict, dfslist, visited, loops, N, originid, -1)
    loops = loops[sortperm(sort.(loops))[1:2:length(loops)]] # removes double entries of loop connections and keeps the first found pair

    return dfsgraph, convert(SVector{N}, dfslist), loops
end

function dfs!(A::Matrix, Adfs::Matrix, dict::Dict, list::Vector, visited::Vector, loops::Vector{<:Vector{<:Integer}},
        index::Integer, currentid::Integer, parentid::Integer
    )

    i = dict[currentid]
    for (childid, j) in dict
        if A[i,j] && parentid != childid # connection from i to j in adjacency && not a direct connection back to the parent
            if visited[j]
                push!(loops, [childid,currentid]) # childid is actually a predecessor of currentid since it's a loop
            else
                index -= 1
                list[index] = childid
                visited[j] = true
                Adfs[i,j] = true
                index = dfs!(A, Adfs, dict, list, visited, loops, index, childid, currentid)
            end
        end
    end
    return index
end

function pattern(dfsgraph::Matrix, dict::Dict, loops::Vector{<:Vector{<:Integer}})
    pat = deepcopy(dfsgraph)

    for loop in loops # loop = [index of starting constraint, index of last body in loop]
        startid = loop[1]
        endid = loop[2]
        currentid = endid
        nextid = parent(dfsgraph, dict, currentid)
        while nextid != startid
            pat[dict[startid],dict[currentid]] = true
            currentid = nextid
            nextid = parent(dfsgraph, dict, nextid)
        end
    end

    return pat
end

function fillins(dfsgraph::Matrix, pattern::Matrix, dict::Dict, loops::Vector{<:Vector{<:Integer}})
    fil = deepcopy(dfsgraph .⊻ pattern) # xor so only fillins remain (+ to be removed loop closure since this is a child)
    originals = deepcopy(dfsgraph)

    for loop in loops # loop = [index of starting constraint, index of last body in loop]
        startid = loop[1]
        endid = loop[2]
        fil[dict[startid],dict[endid]] = false
        originals[dict[startid],dict[endid]] = true
    end

    return convert(Matrix{Bool}, fil), convert(Matrix{Bool}, originals)
end

function parent(dfsgraph::Matrix, dict::Dict, childid::Integer)
    j = dict[childid]
    for (parentid, i) in dict
        dfsgraph[i,j] && (return parentid)
    end
    return -1
end

### Graph functions

# this is done in order!
function successors(dfslist, pattern, dict::Dict)
    N = length(dfslist)
    sucs = [Int64[] for i = 1:N]
    for i = 1:N
        for childid in dfslist
            pattern[i][dict[childid]] && push!(sucs[i], childid)
        end
    end

    return sucs
end

# this is done in order!
function directchildren(dfslist, dfsgraph, dict::Dict)
    N = length(dfslist)
    dirs = [Int64[] for i = 1:N]
    for i = 1:N
        for childid in dfslist
            dfsgraph[i][dict[childid]] && push!(dirs[i], childid)
        end
    end

    return dirs
end

# this is done in order!
function loopchildren(dfslist, fillins, dict::Dict)
    N = length(dfslist)
    loos = [Int64[] for i = 1:N]
    for i = 1:N
        for childid in dfslist
            fillins[i][dict[childid]] && push!(loos[i], childid)
        end
    end

    return loos
end

function ineqchildren(dfslist, bodies, ineqconstraints, dict::Dict)
    N = length(dfslist)
    ineqs = [Int64[] for i = 1:N]
    for body in bodies
        for ineqc in ineqconstraints
            ineqc.parentid == body.id && push!(ineqs[dict[body.id]], ineqc.id)
        end
    end

    return ineqs
end

# this is done in reverse order (but this is not really important for predecessors)
function predecessors(dfslist, pattern, dict::Dict)
    N = length(dfslist)
    preds = [Int64[] for i = 1:N]
    for i = 1:N
        for childid in reverse(dfslist)
            pattern[dict[childid]][i] && push!(preds[i], childid)
        end
    end

    return preds
end

# this is done in order (but this is not really important for connections)
function connections(dfslist, adjacency, dict::Dict)
    N = length(dfslist)
    cons = [Int64[] for i = 1:N]
    for i = 1:N
        for childid in dfslist
            adjacency[i][dict[childid]] && push!(cons[i], childid)
        end
    end

    return cons
end

function recursivedirectchildren!(graph, id::Integer)
    dirs = copy(directchildren(graph, id))
    dirslocal = copy(dirs)
    for childid in dirslocal
        append!(dirs, recursivedirectchildren!(graph, childid))
    end
    return dirs
end

@inline directchildren(graph, id::Integer) = graph.directchildren[graph.dict[id]]
@inline loopchildren(graph, id::Integer) = graph.loopchildren[graph.dict[id]]
@inline ineqchildren(graph, id::Integer) = graph.ineqchildren[graph.dict[id]]
@inline successors(graph, id::Integer) = graph.successors[graph.dict[id]]
@inline predecessors(graph, id::Integer) = graph.predecessors[graph.dict[id]]
@inline connections(graph, id::Integer) = graph.connections[graph.dict[id]]
@inline isactive(graph, id::Integer) = graph.activedict[id]
@inline isinactive(graph, id::Integer) = !isactive(graph, id)

@inline function hassuccessor(graph::Graph, id, childid)
    for val in graph.successors[graph.dict[id]]
        val == childid && (return true)
    end
    return false
end

@inline function haspredecessor(graph::Graph, id, parentid)
    for val in graph.predecessors[graph.dict[id]]
        val == parentid && (return true)
    end
    return false
end

@inline function hasdirectchild(graph::Graph, id, childid)
    for val in graph.directchildren[graph.dict[id]]
        val == childid && (return true)
    end
    return false
end

@inline function activate!(graph::Graph, id::Integer)
    graph.activedict[id] = true
    return
end

@inline function deactivate!(graph::Graph, id::Integer)
    graph.activedict[id] = false
    return
end