How to use the libceed.EVAL_GRAD function in libceed

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

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

    Xq.restore_array_read()
    U.set_array(u, cmode=libceed.USE_POINTER)

    # Calculate G u at quadrature points, G' * 1 at dofs
    bug = ceed.BasisTensorH1Lagrange(dim, 1, P, Q, libceed.GAUSS)
    bug.apply(1, libceed.EVAL_GRAD, U, Uq)
    bug.T.apply(1, libceed.EVAL_GRAD, Ones, Gtposeones)

    # Check if 1' * G * u = u' * (G' * 1)
    gtposeones = Gtposeones.get_array_read()
    uq = Uq.get_array_read()

    for i in range(Pdim):
      sum1 += gtposeones[i]*u[i]
    for i in range(dim*Qdim):
      sum2 += uq[i]
    Gtposeones.restore_array_read()
    Uq.restore_array_read()

    assert math.fabs(sum1 - sum2) < 1E-10

def test_402(ceed_resource, capsys):
  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, "t400-qfunction.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)

  print(qf_setup)
  print(qf_mass)

  stdout, stderr = capsys.readouterr()
  with open(os.path.abspath("./output/test_402.out")) as output_file:
    true_output = output_file.read()

  assert stdout == true_output

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*dim)
  out_vec.set_value(0)

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

  # Check values at quadrature points
  out_array = out_vec.get_array_read()
  for i in range(Q):
    value = dfeval(xq[0*Q+i], xq[1*Q+i])
    assert math.fabs(out_array[0*Q+i] - value) < 1E-10

    value = dfeval(xq[1*Q+i], xq[0*Q+i])
    assert math.fabs(out_array[1*Q+i] - value) < 1E-10

  out_vec.restore_array_read()

qweight = np.empty(q_tet, dtype="float64")
  interp, grad = bm.buildmats(qref, qweight)

  bx_tet = ceed.BasisH1(libceed.TRIANGLE, dim, p_tet, q_hex, interp, grad, qref,
                        qweight)
  bu_tet = ceed.BasisH1(libceed.TRIANGLE, 1, p_tet, q_hex, interp, grad, qref,
                        qweight)

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

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

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

  # Operators
  op_setup_tet = ceed.Operator(qf_setup_tet)
  op_setup_tet.set_field("weights", rxi_tet, bx_tet, libceed.VECTOR_NONE)
  op_setup_tet.set_field("dx", rx_tet, bx_tet, libceed.VECTOR_ACTIVE)
  op_setup_tet.set_field("rho", rui_tet, libceed.BASIS_COLLOCATED,
                         qdata_tet)

  op_mass_tet = ceed.Operator(qf_mass_tet)

indu[p*i+j] = i*(p-1) + j
  ru = ceed.ElemRestriction(nelem, p, nu, 1, indu, cmode=libceed.USE_POINTER)
  rui = ceed.IdentityElemRestriction(nelem, q, q*nelem, 1)

  # Bases
  bx = ceed.BasisTensorH1Lagrange(1, 1, 2, q, libceed.GAUSS)
  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)

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