ENH: Enable querying FEM values from Python

yt/utilities/lib/element_mappings.pyx

@@ -95,6 +95,24 @@ cdef class ElementSampler:
     cdef int check_inside(self, double* mapped_coord) noexcept nogil:
         pass
 
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    @cython.cdivision(True)
+    def check_contains(self, np.float64_t[:,::1] vertices,
+                             np.float64_t[:,::1] positions):
+        cdef np.ndarray[np.float64_t, ndim=2] mapped_coords
+        cdef np.ndarray[np.uint8_t, ndim=1] mask
+        mapped_coords = self.map_reals_to_unit(vertices, positions)
+        mask = np.zeros(mapped_coords.shape[0], dtype=np.uint8) 
+        cdef double[3] mapped_coord
+        cdef int i, j
+        for i in range(mapped_coords.shape[0]):
+            for j in range(mapped_coords.shape[1]):
+                mapped_coord[j] = mapped_coords[i, j]
+            mask[i] = self.check_inside(mapped_coord)
+        return mask
+
+
     @cython.boundscheck(False)
     @cython.wraparound(False)
     @cython.cdivision(True)
@@ -115,6 +133,34 @@ cdef class ElementSampler:
         val = self.sample_at_unit_point(mapped_coord, field_values)
         return val
 
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    @cython.cdivision(True)
+    def map_reals_to_unit(self,
+                          np.float64_t[:,::1] vertices,
+                          np.float64_t[:,::1] positions):
+        cdef double mapped_x[3]
+        cdef int i, n
+        # We have N vertices, which each have three components.
+        cdef np.ndarray[np.float64_t, ndim=2] output_coords
+        output_coords = np.zeros((positions.shape[0], positions.shape[1]), dtype="float64")
+        # Now for each position, we map
+        for n in range(positions.shape[0]):
+            self.map_real_to_unit(mapped_x, &vertices[0,0], &positions[n, 0])
+            for i in range(positions.shape[1]):
+                output_coords[n, i] = mapped_x[i]
+        return output_coords
+
+    def sample_at_real_points(self,
+                              np.float64_t[:,::1] vertices,
+                              np.float64_t[::1] field_values,
+                              np.float64_t[:,::1] positions):
+        cdef np.ndarray[np.float64_t, ndim=1] output_values
+        output_values = np.zeros(positions.shape[0], dtype="float64")
+        for n in range(positions.shape[0]):
+            output_values[n] = self.sample_at_real_point(
+                &vertices[0,0], &field_values[0], &positions[n,0])
+        return output_values
 
 cdef class P1Sampler1D(ElementSampler):
     '''

yt/utilities/mesh_code_generation.py

@@ -83,11 +83,11 @@ def _compute_jacobian(self):
         assert self.num_vertices == len(self.N)
         assert self.num_dim == self.num_mapped_coords
 
-        X = MatrixSymbol("vertices", self.num_vertices, self.num_dim)
+        self.X = MatrixSymbol("vertices", self.num_vertices, self.num_dim)
         self.fx = MatrixSymbol("fx", self.num_dim, 1)
-        physical_position = MatrixSymbol("phys_x", self.num_dim, 1)
+        self.physical_position = MatrixSymbol("phys_x", self.num_dim, 1)
 
-        self.f = (self.N.T * Matrix(X)).T - physical_position
+        self.f = (self.N.T * Matrix(self.X)).T - self.physical_position
 
         self.J = symarray("J", (self.num_dim, self.num_dim))
         for i in range(self.num_dim):