How to use the libceed.EVAL_INTERP function in libceed

for d in range(dim):
      for i in range(Xdim):
        x[d*Xdim + i] = 1 if (i % (2**(dim-d))) // (2**(dim-d-1)) else -1

    X = ceed.Vector(Xdim*dim)
    X.set_array(x, cmode=libceed.USE_POINTER)
    Xq = ceed.Vector(Qdim*dim)
    Xq.set_value(0)
    U = ceed.Vector(Qdim)
    U.set_value(0)
    Uq = ceed.Vector(Qdim)

    bxl = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS_LOBATTO)
    bul = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS_LOBATTO)

    bxl.apply(1, libceed.EVAL_INTERP, X, Xq)

    xq = Xq.get_array_read()
    for i in range(Qdim):
      xx = np.empty(dim, dtype="float64")
      for d in range(dim):
        xx[d] = xq[d*Qdim + i]
      uq[i] = eval(dim, xx)

    Xq.restore_array_read()
    Uq.set_array(uq, cmode=libceed.USE_POINTER)

    # This operation is the identity because the quadrature is collocated
    bul.T.apply(1, libceed.EVAL_INTERP, Uq, U)

    bxg = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS)
    bug = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS)

# Bases
  bx_hex = ceed.BasisTensorH1Lagrange(dim, dim, p_hex, q_hex, libceed.GAUSS)
  bu_hex = ceed.BasisTensorH1Lagrange(dim, 1, p_hex, q_hex, libceed.GAUSS)

  # QFunctions
  qf_setup_hex = ceed.QFunction(1, qfs.setup_mass_2d,
                                os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d"))
  qf_setup_hex.add_input("weights", 1, libceed.EVAL_WEIGHT)
  qf_setup_hex.add_input("dx", dim*dim, libceed.EVAL_GRAD)
  qf_setup_hex.add_output("rho", 1, libceed.EVAL_NONE)

  qf_mass_hex = ceed.QFunction(1, qfs.apply_mass,
                               os.path.join(file_dir, "test-qfunctions.h:apply_mass"))
  qf_mass_hex.add_input("rho", 1, libceed.EVAL_NONE)
  qf_mass_hex.add_input("u", 1, libceed.EVAL_INTERP)
  qf_mass_hex.add_output("v", 1, libceed.EVAL_INTERP)

  # Operators
  op_setup_hex = ceed.Operator(qf_setup_tet)
  op_setup_hex.set_field("weights", rxi_hex, bx_hex, libceed.VECTOR_NONE)
  op_setup_hex.set_field("dx", rx_hex, bx_hex, libceed.VECTOR_ACTIVE)
  op_setup_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED,
                         qdata_hex)

  op_mass_hex = ceed.Operator(qf_mass_hex)
  op_mass_hex.set_field("rho", rui_hex, libceed.BASIS_COLLOCATED, qdata_hex)
  op_mass_hex.set_field("u", ru_hex, bu_hex, libceed.VECTOR_ACTIVE)
  op_mass_hex.set_field("v", ru_hex, bu_hex, libceed.VECTOR_ACTIVE)

  ## ------------------------- Composite Operators -------------------------

interp, grad = bm.buildmats(qref, qweight)

  b = ceed.BasisH1(libceed.TRIANGLE, 1, P, Q, interp, grad, qref, qweight)

  # Interpolate function to quadrature points
  for i in range(P):
    in_array[i] = feval(xr[0*P+i], xr[1*P+i])

  in_vec = ceed.Vector(P)
  in_vec.set_array(in_array, cmode=libceed.USE_POINTER)
  out_vec = ceed.Vector(Q)
  out_vec.set_value(0)
  weights_vec = ceed.Vector(Q)
  weights_vec.set_value(0)

  b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec)
  b.apply(1, libceed.EVAL_WEIGHT, libceed.VECTOR_NONE, weights_vec)

  # Check values at quadrature points
  out_array = out_vec.get_array_read()
  weights_array = weights_vec.get_array_read()
  sum = 0
  for i in range(Q):
    sum += out_array[i]*weights_array[i]
  if math.fabs(sum - 17./24.) > 1E-10:
    print("%f != %f"%(sum, 17./24.))

  out_vec.restore_array_read()
  weights_vec.restore_array_read()

for i in range(Qdim):
      xx = np.empty(dim, dtype="float64")
      for d in range(dim):
        xx[d] = xq[d*Qdim + i]
      uq[i] = eval(dim, xx)

    Xq.restore_array_read()
    Uq.set_array(uq, cmode=libceed.USE_POINTER)

    # This operation is the identity because the quadrature is collocated
    bul.T.apply(1, libceed.EVAL_INTERP, Uq, U)

    bxg = ceed.BasisTensorH1Lagrange(dim, dim, 2, Q, libceed.GAUSS)
    bug = ceed.BasisTensorH1Lagrange(dim, 1, Q, Q, libceed.GAUSS)

    bxg.apply(1, libceed.EVAL_INTERP, X, Xq)
    bug.apply(1, libceed.EVAL_INTERP, U, Uq)

    xq = Xq.get_array_read()
    u = Uq.get_array_read()

    for i in range(Qdim):
      xx = np.empty(dim, dtype="float64")
      for d in range(dim):
        xx[d] = xq[d*Qdim + i]
      fx = eval(dim, xx)
      assert math.fabs(u[i] - fx) < 1E-4

    Xq.restore_array_read()
    Uq.restore_array_read()

def test_412(ceed_resource):
  ceed = libceed.Ceed(ceed_resource)

  size = 3
  qf = ceed.IdentityQFunction(size, libceed.EVAL_INTERP, libceed.EVAL_INTERP)

  q = 8

  u_array = np.zeros(q*size, dtype="float64")
  for i in range(q*size):
    u_array[i] = i*i

  u = ceed.Vector(q*size)
  u.set_array(u_array, cmode=libceed.USE_POINTER)
  v = ceed.Vector(q*size)
  v.set_value(0)

  inputs = [ u ]
  outputs = [ v ]
  qf.apply(q, inputs, outputs)

qweight = np.empty(Q, dtype="float64")

  interp, grad = bm.buildmats(qref, qweight)

  b = ceed.BasisH1(libceed.TRIANGLE, 1, P, Q, interp, grad, qref, qweight)

  # Interpolate function to quadrature points
  for i in range(P):
    in_array[i] = feval(xr[0*P+i], xr[1*P+i])

  in_vec = ceed.Vector(P)
  in_vec.set_array(in_array, cmode=libceed.USE_POINTER)
  out_vec = ceed.Vector(Q)
  out_vec.set_value(0)

  b.apply(1, libceed.EVAL_INTERP, in_vec, out_vec)

  # Check values at quadrature points
  out_array = out_vec.get_array_read()
  for i in range(Q):
    value = feval(xq[0*Q+i], xq[1*Q+i])
    if math.fabs(out_array[i] - value) > 1E-10:
      # LCOV_EXCL_START
      printf("[%d] %f != %f"%(i, out[i], value))
    # LCOV_EXCL_STOP
  out_vec.restore_array_read()

bu = ceed.BasisTensorH1Lagrange(1, 1, p, q, libceed.GAUSS)

  # QFunctions
  file_dir = os.path.dirname(os.path.abspath(__file__))
  qfs = load_qfs_so()

  qf_setup = ceed.QFunction(1, qfs.setup_mass,
                            os.path.join(file_dir, "test-qfunctions.h:setup_mass"))
  qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT)
  qf_setup.add_input("dx", 1, libceed.EVAL_GRAD)
  qf_setup.add_output("rho", 1, libceed.EVAL_NONE)

  qf_mass = ceed.QFunction(1, qfs.apply_mass,
                           os.path.join(file_dir, "test-qfunctions.h:apply_mass"))
  qf_mass.add_input("rho", 1, libceed.EVAL_NONE)
  qf_mass.add_input("u", 1, libceed.EVAL_INTERP)
  qf_mass.add_output("v", 1, libceed.EVAL_INTERP)

  # Operators
  op_setup = ceed.Operator(qf_setup)
  op_setup.set_field("weights", rxi, bx, libceed.VECTOR_NONE)
  op_setup.set_field("dx", rx, bx, libceed.VECTOR_ACTIVE)
  op_setup.set_field("rho", rui, libceed.BASIS_COLLOCATED,
                     libceed.VECTOR_ACTIVE)

  op_mass = ceed.Operator(qf_mass)
  op_mass.set_field("rho", rui, libceed.BASIS_COLLOCATED, qdata)
  op_mass.set_field("u", ru, bu, libceed.VECTOR_ACTIVE)
  op_mass.set_field("v", ru, bu, libceed.VECTOR_ACTIVE)

  # Setup
  op_setup.apply(x, qdata)

def test_400(ceed_resource):
  ceed = libceed.Ceed(ceed_resource)

  file_dir = os.path.dirname(os.path.abspath(__file__))
  qfs = load_qfs_so()

  qf_setup = ceed.QFunction(1, qfs.setup_mass,
                            os.path.join(file_dir, "test-qfunctions.h:setup_mass"))
  qf_setup.add_input("w", 1, libceed.EVAL_WEIGHT)
  qf_setup.add_input("dx", 1, libceed.EVAL_GRAD)
  qf_setup.add_output("qdata", 1, libceed.EVAL_NONE)

  qf_mass = ceed.QFunction(1, qfs.apply_mass,
                           os.path.join(file_dir, "test-qfunctions.h:apply_mass"))
  qf_mass.add_input("qdata", 1, libceed.EVAL_NONE)
  qf_mass.add_input("u", 1, libceed.EVAL_INTERP)
  qf_mass.add_output("v", 1, libceed.EVAL_INTERP)

  q = 8

  w_array = np.zeros(q, dtype="float64")
  u_array = np.zeros(q, dtype="float64")
  v_true  = np.zeros(q, dtype="float64")
  for i in range(q):
    x = 2.*i/(q-1) - 1
    w_array[i] = 1 - x*x
    u_array[i] = 2 + 3*x + 5*x*x
    v_true[i]  = w_array[i] * u_array[i]

  dx = ceed.Vector(q)
  dx.set_value(1)
  w = ceed.Vector(q)