point-in-cell.py

from firedrake import *

mesh = UnitSquareMesh(10, 10)

point_in_cell = """

void point_in_cell(double *out[1], double **coords,
                   double *point)
{
   double x = point[0];
   double y = point[1];

   /* Points of coordinates are

     coords[0][0]; (x0)
     coords[0][1]; (y0)
     ...
     coords[i][0]; (xi)
     coords[i][1]; (yi)

   */

   if ( point is inside cell ) {
        out[0][0] = whatever_you_want;
   }
}
"""

V = FunctionSpace(mesh, 'DG', 0)

point = op2.Global(2, data=[0.5, 0.25])

delta = Function(V)

kernel = op2.Kernel(point_in_cell, "point_in_cell")
op2.par_loop(kernel, delta.cell_set,
             delta.dat(op2.WRITE, delta.cell_node_map()),
             mesh.coordinates.dat(op2.READ, mesh.coordinates.cell_node_map()),
             point(op2.READ))