review fixes to code and documentation

docs/source/point-evaluation.rst

@@ -46,9 +46,15 @@ are equivalent:
    f.at([0.2, 0.4])
    f.at(numpy.array([0.2, 0.4]))
 
-For a single point, the result is ``numpy`` array, or a tuple of
-``numpy`` arrays in case *mixed* functions.  When evaluating multiple
-points, the result is a list of values for each point.
+For a single point, the result is a ``numpy`` array, or a tuple of
+``numpy`` arrays in case of *mixed* functions.  When evaluating
+multiple points, the result is a list of values for each point.
+To summarise:
+
+* Single point, non-mixed: ``numpy`` array
+* Single point, mixed: tuple of ``numpy`` arrays
+* Multiple points, non-mixed: list of ``numpy`` arrays
+* Multiple points, mixed: list of tuples of ``numpy`` arrays
 
 
 Points outside the domain
@@ -57,7 +63,7 @@ Points outside the domain
 When any point is outside the domain of the function,
 :py:class:`.PointNotInDomainError` exception is raised. If
 ``dont_raise=True`` is passed to :func:`~.Function.at`, the result is
-``None``.
+``None`` for those points which fall outside the domain.
 
 .. code-block:: python
 
@@ -73,7 +79,7 @@ When any point is outside the domain of the function,
 
 .. warning::
 
-   Point evaluation on *immersed manifolds* in not supported yet, due
+   Point evaluation on *immersed manifolds* is not supported yet, due
    to the difficulty of specifying a physical point on the manifold.
 
 

firedrake/function.py

@@ -502,7 +502,7 @@ def __idiv__(self, expr):
         return self
 
     @utils.cached_property
-    def ctypes(self):
+    def _ctypes(self):
         # Retrieve data from Python object
         function_space = self.function_space()
         mesh = function_space.mesh()
@@ -533,12 +533,18 @@ def _c_evaluate(self):
 
     def evaluate(self, coord, mapping, component, index_values):
         # Called by UFL when evaluating expressions at coordinates
+        if component or index_values:
+            raise NotImplementedError("Unsupported arguments when attempting to evaluate Function.")
         return self.at(coord)
 
     def at(self, arg, *args, **kwargs):
         """Evaluate function at points."""
         from mpi4py import MPI
-        comm = MPI.COMM_WORLD
+        halo = self.dof_dset.halo
+        if isinstance(halo, tuple):
+            # mixed function space
+            halo = halo[0]
+        comm = halo.comm.tompi4py()
 
         if args:
             arg = (arg,) + args
@@ -573,14 +579,14 @@ def at(self, arg, *args, **kwargs):
 
         def single_eval(x, buf):
             """Helper function to evaluate at a single point."""
-            err = self._c_evaluate(self.ctypes, x.ctypes.data, buf.ctypes.data)
+            err = self._c_evaluate(self._ctypes, x.ctypes.data, buf.ctypes.data)
             if err == -1:
                 raise PointNotInDomainError(self.function_space().mesh(), x.reshape(-1))
 
         if not len(arg.shape) <= 2:
             raise ValueError("Function.at expects point or array of points.")
         points = arg.reshape(-1, arg.shape[-1])
-        value_shape = self.function_space().ufl_element().value_shape()
+        value_shape = self.ufl_shape
 
         split = self.split()
         mixed = len(split) != 1
@@ -609,29 +615,25 @@ def same_result(a, b):
             else:
                 return np.allclose(a, b)
 
-        all_results = comm.gather(l_result, root=0)
-        if comm.rank == 0:
-            g_result = [None] * len(points)
-            for results in all_results:
-                for i, result in results:
-                    if g_result[i] is None:
-                        g_result[i] = result
-                    elif same_result(result, g_result[i]):
-                        pass
-                    else:
-                        raise RuntimeError("Point evaluation gave different results across processes.")
-            result = comm.bcast(g_result, root=0)
-        else:
-            result = comm.bcast(None, root=0)
+        all_results = comm.allgather(l_result)
+        g_result = [None] * len(points)
+        for results in all_results:
+            for i, result in results:
+                if g_result[i] is None:
+                    g_result[i] = result
+                elif same_result(result, g_result[i]):
+                    pass
+                else:
+                    raise RuntimeError("Point evaluation gave different results across processes.")
 
         if not dont_raise:
-            for i in xrange(len(result)):
-                if result[i] is None:
+            for i in xrange(len(g_result)):
+                if g_result[i] is None:
                     raise PointNotInDomainError(self.function_space().mesh(), points[i].reshape(-1))
 
         if len(arg.shape) == 1:
-            result = result[0]
-        return result
+            g_result = g_result[0]
+        return g_result
 
 
 class PointNotInDomainError(Exception):