How to use the libceed.TRIANGLE function in libceed

indx[i*2*p+11] = 16 + offset

  rx = ceed.ElemRestriction(nelem, p, ndofs, dim, indx,
                            cmode=libceed.USE_POINTER)
  rxi = ceed.IdentityElemRestriction(nelem, p, nelem*p, dim)

  ru = ceed.ElemRestriction(nelem, p, ndofs, 1, indx, cmode=libceed.USE_POINTER)
  rui = ceed.IdentityElemRestriction(nelem, q, nqpts, 1)

  # Bases
  qref = np.empty(dim*q, dtype="float64")
  qweight = np.empty(q, dtype="float64")
  interp, grad = bm.buildmats(qref, qweight)

  bx = ceed.BasisH1(libceed.TRIANGLE, dim, p, q, interp, grad, qref, qweight)
  bu = ceed.BasisH1(libceed.TRIANGLE, 1, p, q, interp, grad, qref, qweight)

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

  qf_setup = ceed.QFunction(1, qfs.setup_mass_2d,
                            os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d"))
  qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT)
  qf_setup.add_input("dx", dim*dim, 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)

import sys
import libceed
import numpy as np
import t320buildmats

if __name__ == "__main__":
  ceed = libceed.Ceed(sys.argv[1])

  P, Q, dim = 6, 4, 2

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

  interp, grad = buildmats(qref, qweight)

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

return x1*x1 + x2*x2 + x1*x2 + 1

if __name__ == "__main__":
  ceed = libceed.Ceed(sys.argv[1])

  P, Q, dim = 6, 4, 2

  xq = np.array([0.2, 0.6, 1./3., 0.2, 0.2, 0.2, 1./3., 0.6], dtype="float64")
  xr = np.array([0., 0.5, 1., 0., 0.5, 0., 0., 0., 0., 0.5, 0.5, 1.], dtype="float64")
  in_array = np.empty(P, dtype="float64")
  qref = np.empty(dim*Q, dtype="float64")
  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])

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

  P, Q, dim = 6, 4, 2

  xr = np.array([0., 0.5, 1., 0., 0.5, 0., 0., 0., 0., 0.5, 0.5, 1.], dtype="float64")
  in_array = np.empty(P, dtype="float64")
  qref = np.empty(dim*Q, dtype="float64")
  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)
  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

def feval(x1, x2):
  return x1*x1 + x2*x2 + x1*x2 + 1

if __name__ == "__main__":
  ceed = libceed.Ceed(sys.argv[1])

  P, Q, dim = 6, 4, 2

  xr = np.array([0., 0.5, 1., 0., 0.5, 0., 0., 0., 0., 0.5, 0.5, 1.], dtype="float64")
  in_array = np.empty(P, dtype="float64")
  qref = np.empty(dim*Q, dtype="float64")
  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)
  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

indx[i*2*p+10] =  8 + offset
    indx[i*2*p+11] = 16 + offset

  rx = ceed.ElemRestriction(nelem, p, ndofs, dim, indx,
                            cmode=libceed.USE_POINTER)
  rxi = ceed.IdentityElemRestriction(nelem, p, nelem*p, dim)

  ru = ceed.ElemRestriction(nelem, p, ndofs, 1, indx, cmode=libceed.USE_POINTER)
  rui = ceed.IdentityElemRestriction(nelem, q, nqpts, 1)

  # Bases
  qref = np.empty(dim*q, dtype="float64")
  qweight = np.empty(q, dtype="float64")
  interp, grad = bm.buildmats(qref, qweight)

  bx = ceed.BasisH1(libceed.TRIANGLE, dim, p, q, interp, grad, qref, qweight)
  bu = ceed.BasisH1(libceed.TRIANGLE, 1, p, q, interp, grad, qref, qweight)

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

  qf_setup = ceed.QFunction(1, qfs.setup_mass_2d,
                            os.path.join(file_dir, "test-qfunctions.h:setup_mass_2d"))
  qf_setup.add_input("weights", 1, libceed.EVAL_WEIGHT)
  qf_setup.add_input("dx", dim*dim, 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)