Repository: pyvista/pyvista

Issue:
Feature Request: Glyph filter for cell centered data
**Describe the bug, what's wrong, and what you expect:**

Glyphs can be created for point-based vectors but not cell-based vectors.  It would be nice to support using this with cell based vectors.

-----

**To Reproduce**

This works:

```python
import pyvista as pv
import numpy as np

dims = (10, 10, 10)
spacing = (10, 20, 4)

mesh = pv.UniformGrid(dims, spacing)

n_cells = 9*9*9
n_points = 10*10*10

assert n_cells==mesh.n_cells
assert n_points==mesh.n_points

velocity = np.empty([n_points, 3])
velocity[:] = [0, 1, 0]
mesh["velocity"] = velocity

glyph = mesh.glyph(orient="velocity", factor=2)

print(mesh.active_vectors_info)

plotter = pv.Plotter()
plotter.add_mesh(mesh.outline())
plotter.add_mesh(glyph)

plotter.add_axes()
plotter.show(cpos=[(81.91873552469669, 178.0751383928379, 138.14392189340128),
 (46.0, 89.99999847263098, 18.00000062584877),
 (-0.16676975151033355, 0.8183202262809407, -0.5500362326617142)])
```
![image](https://user-images.githubusercontent.com/39341281/118680711-f5cd0e80-b7cc-11eb-844e-0a3839540234.png)


This doesn't:
```python
import pyvista as pv
import numpy as np

dims = (10, 10, 10)
spacing = (10, 20, 4)

mesh = pv.UniformGrid(dims, spacing)

n_cells = 9*9*9
n_points = 10*10*10

assert n_cells==mesh.n_cells
assert n_points==mesh.n_points

velocity = np.empty([n_cells, 3])
velocity[:] = [0, 1, 0]
mesh["velocity"] = velocity

glyph = mesh.glyph(orient="velocity", factor=2)

print(mesh.active_vectors_info)

plotter = pv.Plotter()
plotter.add_mesh(mesh.outline())
plotter.add_mesh(glyph)

plotter.add_axes()
plotter.show(cpos=[(81.91873552469669, 178.0751383928379, 138.14392189340128),
 (46.0, 89.99999847263098, 18.00000062584877),
 (-0.16676975151033355, 0.8183202262809407, -0.5500362326617142)])

```
![image](https://user-images.githubusercontent.com/39341281/118692939-937a0b00-b7d8-11eb-9a88-98936a371931.png)


----

**System Information:**
 Please run the following code wherever you are experiencing the bug and paste the output below. This report helps us track down bugs and it is critical to addressing your bug:

```py
# Get system info
import pyvista as pv
print(pv.Report())
```

```txt
Python 3.8.2 (default, Jun  2 2020, 16:34:16)  [GCC 4.8.5 20150623 (Red Hat
  4.8.5-36)]

           pyvista : 0.30.1
               vtk : 9.0.1
             numpy : 1.20.1
           imageio : 2.9.0
           appdirs : 1.4.4
            scooby : 0.5.6
            meshio : 4.3.9
        matplotlib : 3.3.4
             PyQt5 : 5.15.4
           IPython : 7.21.0
             scipy : 1.6.0
              tqdm : 4.58.0

```


Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/pyvista/pyvista/pull/1334?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pyvista) Report
> Merging [#1334](https://codecov.io/gh/pyvista/pyvista/pull/1334?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pyvista) (6ae99c7) into [master](https://codecov.io/gh/pyvista/pyvista/commit/d6ada7bbd0db3e03605a8736e27dd70b95dddcb5?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pyvista) (d6ada7b) will **decrease** coverage by `0.06%`.
> The diff coverage is `100.00%`.

```diff
@@            Coverage Diff             @@
##           master    #1334      +/-   ##
==========================================
- Coverage   90.05%   89.98%   -0.07%     
==========================================
  Files          56       56              
  Lines       12181    12190       +9     
==========================================
  Hits        10969    10969              
- Misses       1212     1221       +9     
```

MatthewFlamm: The falling test passes locally, so I need to do some investigating.
adeak: Do you have access to the CI results? I also can't reproduce the failure locally, even though I'm on linux and both linux and mac CI fails.
MatthewFlamm: I only tested this one year at a time, so I suspect that one of the other tests is modifying the datasets.  If true, we should move this to a fixture.
MatthewFlamm: > Do you have access to the CI results?

I do, but the failure isn't very helpful.  I'm going to make the dataset creation a fixture, so that it better separates the tests first.


MatthewFlamm: I _think_ it is fixed with mixing to a fixture, as I ran it in my fork.   But this `DATASET` variable is imported to other tests. i will have to go in and clean up there too.
MatthewFlamm: When the latest commit ran tests in my fork, it failed this test, which is untouched here, although I did modify a test in that module.

```txt
FAILED tests/test_common.py::test_rotate_should_match_vtk_rotation[z]
```

It passed in this CI.  

Patch:
diff --git a/docs/extras/building_vtk.rst b/docs/extras/building_vtk.rst
index f466e611..cce56b9f 100644
--- a/docs/extras/building_vtk.rst
+++ b/docs/extras/building_vtk.rst
@@ -2,78 +2,256 @@
 
 Building VTK
 ============
-Kitware provides Python wheels for VTK at `PyPI VTK <https://pypi.org/project/vtk/>`_, but there are situations where you
+Kitware provides Python wheels for VTK at `PyPI VTK
+<https://pypi.org/project/vtk/>`_, but there are situations where you
 may need to build VTK from source (e.g. new release of Python, EGL
 rendering, additional features, etc).  As ``pyvista`` does not provide
 ``vtk``, you will have to build it manually.
 
+Should you need a prebuilt wheel, a variety of prebuilt wheels
+(including ones with EGL) can be found at `pyvista Release 0.31.0
+<https://github.com/pyvista/pyvista/releases/tag/0.31.0>`_.  This
+includes wheels for Raspberry Pi (64-bit OS), Python 3.9 Linux, and
+EGL support (``arch64`` and ``x86_64``).
+
 
 Wheels on Linux and Mac OS
 ~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-Building VTK from source on Linux is fairly easy.  Using the default
-build settings, build a Python wheel of VTK using ``ninja`` using the following script.  This script assumes Python 3.9, but you can use any modern Python version.
-For some additional useful options, see the `conda-forge recipe <https://github.com/conda-forge/vtk-feedstock/blob/master/recipe/build.sh>`__.
-Most of the ones below are designed to reduce the build time and resulting wheel size.
+Building VTK from source on Linux and MacOS is fairly easy.  Using the
+default build settings, build a Python wheel of VTK using ``ninja``
+using the following script.  This script uses system python3, but you
+can use any modern Python version.  For some additional useful
+options, see the `conda-forge recipe
+<https://github.com/conda-forge/vtk-feedstock/blob/master/recipe/build.sh>`__.
+Most of the ones below are designed to reduce the build time and
+resulting wheel size.
 
 .. code-block:: bash
 
+    #!/bin/bash
+
+    # install build dependencies (Linux/Debian)
+    sudo apt-get install ninja-build cmake libgl1-mesa-dev python3-dev -y
+
+    # Linux/CentOS
+    sudo yum install epel-release -y
+    sudo yum install ninja-build cmake mesa-libGL-devel mesa-libGLU-devel -y
+
     git clone https://github.com/Kitware/VTK
     cd VTK
 
-    mkdir build
-    cd build
-    PYBIN=/usr/bin/python3.9
-    cmake -GNinja -DCMAKE_BUILD_TYPE=Release -DVTK_BUILD_TESTING=OFF \
-        -DVTK_BUILD_DOCUMENTATION=OFF -DVTK_BUILD_EXAMPLES=OFF -DVTK_DATA_EXCLUDE_FROM_ALL:BOOL=ON \
-        -DVTK_MODULE_ENABLE_VTK_PythonInterpreter:STRING=NO \
-        -DVTK_WHEEL_BUILD=ON -DVTK_PYTHON_VERSION=3 -DVTK_WRAP_PYTHON=ON \
-        -DPython3_EXECUTABLE=$PYBIN ../
+    PYBIN=/usr/bin/python3
+    cmake -GNinja \
+          -DCMAKE_BUILD_TYPE=Release \
+          -DVTK_BUILD_TESTING=OFF \
+          -DVTK_BUILD_DOCUMENTATION=OFF \
+          -DVTK_BUILD_EXAMPLES=OFF \
+          -DVTK_DATA_EXCLUDE_FROM_ALL:BOOL=ON \
+          -DVTK_MODULE_ENABLE_VTK_PythonInterpreter:STRING=NO \
+          -DVTK_WHEEL_BUILD=ON \
+          -DVTK_PYTHON_VERSION=3 \
+          -DVTK_WRAP_PYTHON=ON \
+          -DVTK_OPENGL_HAS_EGL=False \
+          -DPython3_EXECUTABLE=$PYBIN ../
     ninja
     $PYBIN setup.py bdist_wheel
     pip install dist/vtk-*.whl  # optionally install it
 
-You may need to install ``python3.9-dev`` and ``ninja`` if you have
-not already installed it.
-
-
 .. _gpu_off_screen:
 
+
 Off-Screen Plotting GPU Support
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+VTK supports rendering with EGL, enabling rapid off-screen using GPU
+hardware acceleration.  The default VTK wheels are not built with this
+feature, but you can build VTK for off-screen plotting using GPU
+support by modifying the above ``cmake`` command with:
 
-VTK supports rendering with EGL, enabling rapid off-screen
-using GPU hardware acceleration.  Unfortuantly, the default VTK wheels
-are not built with this feature as it results in a > 400 MB wheel.
-For the adventurous/desperate, build VTK with EGL for a given Python wheel on
-Linux with the following:
+.. code-block:: bash
 
-You can build VTK for off-screen plotting using GPU support by modifying
-the above ``cmake`` command with:
+    # install OpenGL with EGL (Linux/Debian)
+    sudo apt-get install libegl1-mesa-dev -y
 
-.. code-block:: bash
+    # CentOS
+    sudo yum install mesa-libEGL-devel -y
 
-  cmake -GNinja \
-    -DVTK_BUILD_TESTING=OFF \
-    -DVTK_WHEEL_BUILD=ON \
-    -DVTK_PYTHON_VERSION=3 \
-    -DVTK_WRAP_PYTHON=ON \
     -DVTK_OPENGL_HAS_EGL=True \
-    -DVTK_USE_X=False \
-    -DPython3_EXECUTABLE=$PYBIN \
-    ../
-  ninja
-
-  $PYBIN setup.py bdist_wheel
-  $PYBIN -m pip install dist/vtk-*.whl  # optional
 
 This disables any plotting using the X server, so be prepared to use
 this module only on a headless display where you either intend to save
 static images or stream the render window to another computer with a
-display (e.g using ``use_ipyvtk=True`` and jupyterlab). In other words,
-this wheel will make VTK unusable outside of an off-screen
-environment, so only plan on installing it on a headless system
-without a X server.
+display (e.g using ``pyvista.set_jupyter_backend('ipyvtklink')`` and
+jupyterlab). In other words, this wheel will make VTK unusable outside
+of an off-screen environment, so only plan on installing it on a
+headless system without a X server.
+
+
+Building ManyLinux Wheels
+~~~~~~~~~~~~~~~~~~~~~~~~~
+The above directions are great for building a local build of VTK, but
+these wheels are difficult to share outside your local install given
+issues with ABI compatibility due to the version of Linux they were
+built on.  You can work around this by building your wheels using a
+`manylinux <https://github.com/pypa/manylinux>`_ docker image.
+
+.. code-block:: bash
+
+    #!/bin/bash
+    # builds python wheels on docker container and tests installation
+
+    set -e -x
+
+    # build based on python version from args
+    PYTHON_VERSION="$1"
+    case $PYTHON_VERSION in
+    2.7)
+      PYBIN="/opt/python/cp27-cp27m/bin/python"
+      ;;
+    3.5)
+      PYBIN="/opt/python/cp35-cp35m/bin/python"
+      ;;
+    3.6)
+      PYBIN="/opt/python/cp36-cp36m/bin/python"
+      ;;
+    3.7)
+      PYBIN="/opt/python/cp37-cp37m/bin/python"
+      ;;
+    3.8)
+      PYBIN="/opt/python/cp38-cp38/bin/python"
+      ;;
+    3.9)
+      PYBIN="/opt/python/cp39-cp39/bin/python"
+      ;;
+    esac
+
+    yum install ninja-build cmake mesa-libGL-devel mesa-libGLU-devel -y
+
+    rm -rf /io/build
+    mkdir /io/build -p
+    cd /io/build
+
+    cmake -GNinja \
+          -DCMAKE_BUILD_TYPE=Release \
+          -DVTK_BUILD_TESTING=OFF \
+          -DVTK_BUILD_DOCUMENTATION=OFF \
+          -DVTK_BUILD_EXAMPLES=OFF \
+          -DVTK_DATA_EXCLUDE_FROM_ALL:BOOL=ON \
+          -DVTK_MODULE_ENABLE_VTK_PythonInterpreter:STRING=NO \
+          -DVTK_WHEEL_BUILD=ON \
+          -DVTK_PYTHON_VERSION=3 \
+          -DVTK_WRAP_PYTHON=ON \
+          -DVTK_OPENGL_HAS_EGL=False \
+          -DPython3_EXECUTABLE=$PYBIN ../
+    ninja-build
+
+    # build wheel
+    rm -rf dist
+    $PYBIN setup.py bdist_wheel
+
+    # cleanup wheel
+    rm -rf wheelhouse
+    auditwheel repair dist/*.whl
+    cp wheelhouse/vtk*.whl /io/wheels
+
+This script can be called with:
+
+.. code::
+
+    PYTHON_VERSION=3.8
+    docker run -e \
+           --rm -v `pwd`:/io quay.io/pypa/manylinux2014_x86_64 \
+           /io/build_wheels.sh $PYTHON_VERSION
+
+.. note::
+   To build the EGL version of the wheel, follow the directions in the
+   previous section.  Add ``mesa-libEGL-devel`` to the installation
+   dependencies.
+
+
+Python VTK Wheel on Raspberry PI (64-bit)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+While it's possible to build on 32-bit Raspberry Pi (ARMv7), there's
+several issues that crop up when building wheels for the 32-bit
+version (see `manylinux issue 84
+<https://github.com/pypa/manylinux/issues/84>`_.  Should you attempt
+to build on 32-bit , attempt to build the wheel using `dockcross
+<https://github.com/dockcross/dockcross>`_ as you may run into memory limitations otherwise (especially with only 1 GB RAM).
+
+Building the ``aarch64`` manylinux wheel can be done via docker with the ``quay.io/pypa/manylinux2014_aarch64`` image.  Run the following:
+
+.. code::
+
+    PYTHON_VERSION=3.6
+    rm -rf build
+    docker run -e \
+           --rm -v `pwd`:/io quay.io/pypa/manylinux2014_aarch64 \
+           /io/build_wheels.sh $PYTHON_VERSION
+
+Where ``build_wheels.sh`` is:
+
+.. code::
+
+
+    #!/bin/bash
+    # builds python wheels on docker container and tests installation
+
+    set -e -x
+
+    # build based on python version from args
+    PYTHON_VERSION="$1"
+    case $PYTHON_VERSION in
+    2.7)
+      PYBIN="/opt/python/cp27-cp27m/bin/python"
+      ;;
+    3.5)
+      PYBIN="/opt/python/cp35-cp35m/bin/python"
+      ;;
+    3.6)
+      PYBIN="/opt/python/cp36-cp36m/bin/python"
+      ;;
+    3.7)
+      PYBIN="/opt/python/cp37-cp37m/bin/python"
+      ;;
+    3.8)
+      PYBIN="/opt/python/cp38-cp38/bin/python"
+      ;;
+    3.9)
+      PYBIN="/opt/python/cp39-cp39/bin/python"
+      ;;
+    esac
+
+    /bin/bash
+    yum install epel-release -y
+    yum install ninja-build mesa-libEGL-devel -y
+
+    mkdir /io/build -p
+    cd /io/build
+
+    cmake -GNinja \
+          -DCMAKE_BUILD_TYPE=Release \
+          -DVTK_BUILD_TESTING=OFF \
+          -DVTK_BUILD_DOCUMENTATION=OFF \
+          -DVTK_BUILD_EXAMPLES=OFF \
+          -DVTK_DATA_EXCLUDE_FROM_ALL:BOOL=ON \
+          -DVTK_MODULE_ENABLE_VTK_PythonInterpreter:STRING=NO \
+          -DVTK_WHEEL_BUILD=ON \
+          -DVTK_PYTHON_VERSION=3 \
+          -DVTK_WRAP_PYTHON=ON \
+          -DVTK_OPENGL_HAS_EGL=True \
+          -DPython3_EXECUTABLE=$PYBIN ../
+    ninja-build
+
+    # build wheel
+    rm -rf dist
+    $PYBIN setup.py bdist_wheel
+
+    # cleanup wheel
+    rm -rf wheelhouse
+    auditwheel repair dist/*.whl
+    cp wheelhouse/vtk*.whl /io/wheels_egl
+
+Be sure to either enable or disable ``DVTK_OPENGL_HAS_EGL`` depending
+on if you want ``EGL`` enabled for your wheel.
 
 
 Building VTK on Windows
diff --git a/docs/plotting/plotting.rst b/docs/plotting/plotting.rst
index 894e4c8a..2a938325 100644
--- a/docs/plotting/plotting.rst
+++ b/docs/plotting/plotting.rst
@@ -13,6 +13,8 @@ shortcuts are available:
 +-------------------------------------+-----------------+-----------------------------------------------------+
 | ``q``                                                 | Close the rendering window                          |
 +-------------------------------------+-----------------+-----------------------------------------------------+
+| ``f``                                                 | Focus and zoom in on a point                        |
++-------------------------------------+-----------------+-----------------------------------------------------+
 | ``v``                                                 | Isometric camera view                               |
 +-------------------------------------+-----------------+-----------------------------------------------------+
 | ``w``                                                 | Switch all datasets to a `wireframe` representation |
diff --git a/pyvista/core/filters.py b/pyvista/core/filters.py
index f3dab7bb..d004e8a9 100644
--- a/pyvista/core/filters.py
+++ b/pyvista/core/filters.py
@@ -331,6 +331,7 @@ class DataSetFilters:
         >>> from pyvista import examples
         >>> dataset = examples.load_hexbeam()
         >>> clipped = dataset.clip_scalar(value=100, invert=False)
+
         """
         if isinstance(dataset, _vtk.vtkPolyData):
             alg = _vtk.vtkClipPolyData()
@@ -633,6 +634,7 @@ class DataSetFilters:
         >>> volume[:3] = 1
         >>> v = pyvista.wrap(volume)
         >>> threshed = v.threshold(0.1)
+
         """
         # set the scalaras to threshold on
         if scalars is None:
@@ -1068,10 +1070,12 @@ class DataSetFilters:
         Map a puppy texture to a sphere
 
         >>> import pyvista
+        >>> from pyvista import examples
         >>> sphere = pyvista.Sphere()
         >>> sphere = sphere.texture_map_to_sphere()
         >>> tex = examples.download_puppy_texture()  # doctest:+SKIP
         >>> cpos = sphere.plot(texture=tex)  # doctest:+SKIP
+
         """
         alg = _vtk.vtkTextureMapToSphere()
         if center is None:
@@ -1143,23 +1147,27 @@ class DataSetFilters:
     def glyph(dataset, orient=True, scale=True, factor=1.0, geom=None,
               indices=None, tolerance=None, absolute=False, clamping=False,
               rng=None, progress_bar=False):
-        """Copy a geometric representation (called a glyph) to every point in the input dataset.
+        """Copy a geometric representation (called a glyph) to the input dataset.
 
         The glyph may be oriented along the input vectors, and it may
         be scaled according to scalar data or vector
         magnitude. Passing a table of glyphs to choose from based on
-        scalars or vector magnitudes is also supported.
+        scalars or vector magnitudes is also supported.  The arrays
+        used for ``orient`` and ``scale`` must be either both point data
+        or both cell data.
 
         Parameters
         ----------
-        orient : bool
-            Use the active vectors array to orient the glyphs
+        orient : bool or str, optional
+            If ``True``, use the active vectors array to orient the glyphs.
+            If string, the vector array to use to orient the glyphs.
 
-        scale : bool
-            Use the active scalars to scale the glyphs
+        scale : bool or str, optional
+            If ``True``, use the active scalars to scale the glyphs.
+            If string, the scalar array to use to scale the glyphs.
 
-        factor : float
-            Scale factor applied to scaling array
+        factor : float, optional
+            Scale factor applied to scaling array.
 
         geom : vtk.vtkDataSet or tuple(vtk.vtkDataSet), optional
             The geometry to use for the glyph. If missing, an arrow glyph
@@ -1180,13 +1188,13 @@ class DataSetFilters:
             Specify tolerance in terms of fraction of bounding box length.
             Float value is between 0 and 1. Default is None. If ``absolute``
             is ``True`` then the tolerance can be an absolute distance.
-            If None, points merging as a preprocessing step is disabled.
+            If ``None``, points merging as a preprocessing step is disabled.
 
         absolute : bool, optional
             Control if ``tolerance`` is an absolute distance or a fraction.
 
-        clamping: bool
-            Turn on/off clamping of "scalar" values to range.
+        clamping: bool, optional
+            Turn on/off clamping of "scalar" values to range. Default ``False``.
 
         rng: tuple(float), optional
             Set the range of values to be considered by the filter when scalars
@@ -1195,6 +1203,21 @@ class DataSetFilters:
         progress_bar : bool, optional
             Display a progress bar to indicate progress.
 
+        Returns
+        -------
+        glyphs : pyvista.PolyData
+            Glyphs at either the cell centers or points.
+
+        Examples
+        --------
+        Create arrow glyphs oriented by vectors and scaled by scalars.
+        Factor parameter is used to reduce the size of the arrows.
+
+        >>> import pyvista as pv
+        >>> from pyvista import examples
+        >>> mesh = examples.download_carotid().threshold(145, scalars="scalars")  # doctest:+SKIP
+        >>> glyph = mesh.glyph(orient="vectors", scale="scalars", factor=0.01)  # doctest:+SKIP
+
         """
         # Clean the points before glyphing
         if tolerance is not None:
@@ -1253,10 +1276,24 @@ class DataSetFilters:
         if isinstance(orient, str):
             dataset.active_vectors_name = orient
             orient = True
+        if scale and orient:
+            if (dataset.active_vectors_info.association == FieldAssociation.CELL
+                and dataset.active_scalars_info.association == FieldAssociation.CELL
+            ):
+                source_data = dataset.cell_centers()
+            elif(dataset.active_vectors_info.association == FieldAssociation.POINT
+                and dataset.active_scalars_info.association == FieldAssociation.POINT
+            ):
+                source_data = dataset
+            else:
+                raise ValueError("Both ``scale`` and ``orient`` must use "
+                                 "point data or cell data.")
+        else:
+            source_data = dataset
         if rng is not None:
             alg.SetRange(rng)
         alg.SetOrient(orient)
-        alg.SetInputData(dataset)
+        alg.SetInputData(source_data)
         alg.SetVectorModeToUseVector()
         alg.SetScaleFactor(factor)
         alg.SetClamping(clamping)
@@ -1566,6 +1603,7 @@ class DataSetFilters:
 
         progress_bar : bool, optional
             Display a progress bar to indicate progress.
+
         """
         alg = _vtk.vtkDelaunay3D()
         alg.SetInputData(dataset)
@@ -1770,9 +1808,9 @@ class DataSetFilters:
 
         n_points : int, optional
             If given, specifies the number of the closest points used to form
-            the interpolation basis. This will invalidate the radius and
-            sharpness arguments in favor of an N closest points approach. This
-            typically has poorer results.
+            the interpolation basis. This will invalidate the radius argument
+            in favor of an N closest points approach. This typically has poorer
+            results.
 
         strategy : str, optional
             Specify a strategy to use when encountering a "null" point during
@@ -1896,6 +1934,7 @@ class DataSetFilters:
         source : pyvista.PolyData
             The points of the source are the seed points for the streamlines.
             Only returned if ``return_source=True``.
+
         """
         if source_center is None:
             source_center = dataset.center
@@ -1938,8 +1977,7 @@ class DataSetFilters:
                     max_steps=2000, terminal_speed=1e-12, max_error=1e-6,
                     max_time=None, compute_vorticity=True, rotation_scale=1.0,
                     interpolator_type='point'):
-        """
-        Generate streamlines of vectors from the points of a source mesh.
+        """Generate streamlines of vectors from the points of a source mesh.
         
         The integration is performed using a specified integrator, by default
         Runge-Kutta2. This supports integration through any type of dataset.
@@ -1947,8 +1985,8 @@ class DataSetFilters:
         ``surface_streamlines`` parameter is used, the integration is constrained
         to lie on the surface defined by 2D cells.
 
-        Parameters:
-        -----------
+        Parameters
+        ----------
         source : pyvista.DataSet
             The points of the source provide the starting points of the
             streamlines.  This will override both sphere and line sources.
@@ -2024,6 +2062,7 @@ class DataSetFilters:
             (i.e., polyline) representing a streamline. The attribute values
             associated with each streamline are stored in the cell data, whereas
             those associated with streamline-points are stored in the point data.
+
         """
         integration_direction = str(integration_direction).strip().lower()
         if integration_direction not in ['both', 'back', 'backward', 'forward']:
@@ -2126,6 +2165,7 @@ class DataSetFilters:
         -------
         sampled_line : pv.PolyData
             Line object with sampled data from dataset.
+
         """
         if resolution is None:
             resolution = int(dataset.n_cells)
@@ -2258,6 +2298,7 @@ class DataSetFilters:
         >>> pointb = [uniform.bounds[1], uniform.bounds[2], uniform.bounds[4]]
         >>> center = [uniform.bounds[0], uniform.bounds[2], uniform.bounds[4]]
         >>> sampled_arc = uniform.sample_over_circular_arc(pointa, pointb, center)
+
         """
         if resolution is None:
             resolution = int(dataset.n_cells)
@@ -2312,6 +2353,7 @@ class DataSetFilters:
         >>> polar = [uniform.bounds[0], uniform.bounds[2], uniform.bounds[5]]
         >>> center = [uniform.bounds[0], uniform.bounds[2], uniform.bounds[4]]
         >>> sampled_arc = uniform.sample_over_circular_arc_normal(center, normal=normal, polar=polar)
+
         """
         if resolution is None:
             resolution = int(dataset.n_cells)
@@ -2388,6 +2430,7 @@ class DataSetFilters:
         >>> b = [mesh.bounds[1], mesh.bounds[2], mesh.bounds[4]]
         >>> center = [mesh.bounds[0], mesh.bounds[2], mesh.bounds[4]]
         >>> mesh.plot_over_circular_arc(a, b, center, resolution=1000, show=False)
+
         """
         # Ensure matplotlib is available
         try:
@@ -3045,6 +3088,7 @@ class DataSetFilters:
         >>> import pyvista
         >>> mesh = pyvista.Sphere()
         >>> shrunk_mesh = mesh.shrink(shrink_factor=0.8)
+
         """
         if not (0.0 <= shrink_factor <= 1.0):
             raise ValueError('`shrink_factor` should be between 0.0 and 1.0')
@@ -3092,6 +3136,7 @@ class DataSetFilters:
         ...                              [0, 0, 0, 1]])
         >>> transformed = mesh.transform(transform_matrix)
         >>> cpos = transformed.plot(show_edges=True)
+
         """
         if isinstance(trans, _vtk.vtkMatrix4x4):
             m = trans
@@ -3592,6 +3637,7 @@ class PolyDataFilters(DataSetFilters):
         >>> from pyvista import examples
         >>> hills = examples.load_random_hills()
         >>> cpos = hills.plot_curvature(smooth_shading=True)
+
         """
         kwargs.setdefault('scalar_bar_args',
                           {'title': f'{curv_type.capitalize()} Curvature'})
@@ -3682,6 +3728,7 @@ class PolyDataFilters(DataSetFilters):
         Sharp Edges on Cube:        384
         >>> print(f'Sharp Edges on Smooth Cube: {n_smooth_cells}')
         Sharp Edges on Smooth Cube: 12
+
         """
         alg = _vtk.vtkSmoothPolyDataFilter()
         alg.SetInputData(poly_data)
@@ -4614,6 +4661,7 @@ class PolyDataFilters(DataSetFilters):
         ... string = ", ".join([f"({point[0]:.3f}, {point[1]:.3f}, {point[2]:.3f})" for point in points])
         ... print(f'Rays intersected at {string}')
         Rays intersected at (0.499, 0.000, 0.000), (0.000, 0.497, 0.000), (0.000, 0.000, 0.500)
+
         """
         if not poly_data.is_all_triangles():
             raise NotAllTrianglesError
@@ -4626,29 +4674,44 @@ class PolyDataFilters(DataSetFilters):
                 "\tconda install trimesh rtree pyembree"
             )
 
+        origins = np.asarray(origins)
+        directions = np.asarray(directions)
         faces_as_array = poly_data.faces.reshape((poly_data.n_faces, 4))[:, 1:]
         tmesh = trimesh.Trimesh(poly_data.points, faces_as_array)
         locations, index_ray, index_tri = tmesh.ray.intersects_location(
             origins, directions, multiple_hits=not first_point
         )
         if retry:
-            ray_tuples = [(id_r, l, id_t) for id_r, l, id_t in zip(index_ray, locations, index_tri)]
-            for id_r in range(len(origins)):
-                if id_r not in index_ray:
-                    origin = np.array(origins[id_r])
-                    vector = np.array(directions[id_r])
-                    unit_vector = vector / np.sqrt(np.sum(np.power(vector, 2)))
-                    second_point = origin + (unit_vector * poly_data.length)
-                    locs, indexes = poly_data.ray_trace(origin, second_point, first_point=first_point)
-                    if locs.any():
-                        if first_point:
-                            locs = locs.reshape([1, 3])
-                        for loc, id_t in zip(locs, indexes):
-                            ray_tuples.append((id_r, loc, id_t))
-            sorted_results = sorted(ray_tuples, key=lambda x: x[0])
-            locations = np.array([loc for id_r, loc, id_t in sorted_results])
-            index_ray = np.array([id_r for id_r, loc, id_t in sorted_results])
-            index_tri = np.array([id_t for id_r, loc, id_t in sorted_results])
+            # gather intersecting rays in lists
+            loc_lst, ray_lst, tri_lst = [arr.tolist() for arr in [locations, index_ray, index_tri]]
+
+            # find indices that trimesh failed on
+            all_ray_indices = np.arange(len(origins))
+            retry_ray_indices = np.setdiff1d(all_ray_indices, index_ray, assume_unique=True)
+
+            # compute ray points for all failed rays at once
+            origins_retry = origins[retry_ray_indices, :]  # shape (n_retry, 3)
+            directions_retry = directions[retry_ray_indices, :]
+            unit_directions = directions_retry / np.linalg.norm(directions_retry,
+                                                                axis=1, keepdims=True)
+            second_points = origins_retry + unit_directions * poly_data.length  # shape (n_retry, 3)
+
+            for id_r, origin, second_point in zip(retry_ray_indices, origins_retry, second_points):
+                locs, indices = poly_data.ray_trace(origin, second_point, first_point=first_point)
+                if locs.any():
+                    if first_point:
+                        locs = locs.reshape([1, 3])
+                    ray_lst.extend([id_r] * indices.size)
+                    tri_lst.extend(indices)
+                    loc_lst.extend(locs)
+
+            # sort result arrays by ray index
+            index_ray = np.array(ray_lst)
+            sorting_inds = index_ray.argsort()
+            index_ray = index_ray[sorting_inds]
+            index_tri = np.array(tri_lst)[sorting_inds]
+            locations = np.array(loc_lst)[sorting_inds]
+
         return locations, index_ray, index_tri
 
     def plot_boundaries(poly_data, edge_color="red", **kwargs):
@@ -5118,6 +5181,7 @@ class PolyDataFilters(DataSetFilters):
         >>> arc = pyvista.CircularArc([-1, 0, 0], [1, 0, 0], [0, 0, 0])
         >>> mesh = arc.extrude([0, 0, 1])
         >>> cpos = mesh.plot()
+
         """
         alg = _vtk.vtkLinearExtrusionFilter()
         alg.SetExtrusionTypeToVectorExtrusion()
@@ -5191,6 +5255,7 @@ class PolyDataFilters(DataSetFilters):
         >>> line = pyvista.Line(pointa=(0, 0, 0), pointb=(1, 0, 0))
         >>> mesh = line.extrude_rotate(resolution = 4)
         >>> cpos = mesh.plot()
+
         """
         if resolution <= 0:
             raise ValueError('`resolution` should be positive')
@@ -5265,6 +5330,7 @@ class PolyDataFilters(DataSetFilters):
         >>> stripped = slc.strip()
         >>> stripped.n_cells
         1
+
         """
         alg = _vtk.vtkStripper()
         alg.SetInputDataObject(poly_data)
@@ -5292,6 +5358,7 @@ class UnstructuredGridFilters(DataSetFilters):
         ----------
         progress_bar : bool, optional
             Display a progress bar to indicate progress.
+
         """
         return pyvista.PolyData(ugrid.points).delaunay_2d(tol=tol, alpha=alpha,
                                                           offset=offset,
@@ -5352,6 +5419,7 @@ class StructuredGridFilters(DataSetFilters):
 
         >>> joined = voi_1.concatenate(voi_2, axis=1)
         >>> assert np.allclose(grid.points, joined.points)
+
         """
         alg = _vtk.vtkExtractGrid()
         alg.SetVOI(voi)
@@ -5398,6 +5466,7 @@ class StructuredGridFilters(DataSetFilters):
         >>> joined = voi_1.concatenate(voi_2, axis=1)
         >>> print(grid.dimensions, 'same as', joined.dimensions)
         [80, 80, 1] same as [80, 80, 1]
+
         """
         if axis > 2:
             raise RuntimeError('Concatenation axis must be <= 2.')
@@ -5495,6 +5564,7 @@ class UniformGridFilters(DataSetFilters):
 
         progress_bar : bool, optional
             Display a progress bar to indicate progress.
+
         """
         alg = _vtk.vtkImageGaussianSmooth()
         alg.SetInputDataObject(dataset)
@@ -5544,6 +5614,7 @@ class UniformGridFilters(DataSetFilters):
             this is on, the subsampling will always include the boundary of
             the grid even though the sample rate is not an even multiple of
             the grid dimensions. (By default this is off.)
+
         """
         alg = _vtk.vtkExtractVOI()
         alg.SetVOI(voi)


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index 4267f011..b128024f 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -76,3 +76,13 @@ def spline():
 def tri_cylinder():
     """Triangulated cylinder"""
     return pyvista.Cylinder().triangulate()
+
+@fixture()
+def datasets():
+    return [
+        examples.load_uniform(),  # UniformGrid
+        examples.load_rectilinear(),  # RectilinearGrid
+        examples.load_hexbeam(),  # UnstructuredGrid
+        examples.load_airplane(),  # PolyData
+        examples.load_structured(),  # StructuredGrid
+    ]
diff --git a/tests/test_common.py b/tests/test_common.py
index c5e864fa..9a709cee 100644
--- a/tests/test_common.py
+++ b/tests/test_common.py
@@ -7,8 +7,6 @@ from hypothesis.extra.numpy import arrays, array_shapes
 from hypothesis.strategies import composite, integers, floats, one_of
 from vtk.util.numpy_support import vtk_to_numpy
 
-from .test_filters import DATASETS
-
 import pyvista
 from pyvista import examples, Texture
 
@@ -992,32 +990,37 @@ def test_cell_type(grid):
     assert isinstance(ctype, int)
 
 
-@pytest.mark.parametrize('dataset', DATASETS)
-def test_serialize_deserialize(dataset):
-    dataset_2 = pickle.loads(pickle.dumps(dataset))
 
-    # check python attributes are the same
-    for attr in dataset.__dict__:
-        assert getattr(dataset_2, attr) == getattr(dataset, attr)
+def test_serialize_deserialize(datasets):
+    for dataset in datasets:
+        dataset_2 = pickle.loads(pickle.dumps(dataset))
 
-    # check data is the same
-    for attr in ('n_cells', 'n_points', 'n_arrays'):
-        if hasattr(dataset, attr):
+        # check python attributes are the same
+        for attr in dataset.__dict__:
             assert getattr(dataset_2, attr) == getattr(dataset, attr)
 
-    for attr in ('cells', 'points'):
-        if hasattr(dataset, attr):
-            assert getattr(dataset_2, attr) == \
-                   pytest.approx(getattr(dataset, attr))
-
-    for name in dataset.point_arrays:
-        assert dataset_2.point_arrays[name] == \
-               pytest.approx(dataset.point_arrays[name])
-
-    for name in dataset.cell_arrays:
-        assert dataset_2.cell_arrays[name] == \
-               pytest.approx(dataset.cell_arrays[name])
-
-    for name in dataset.field_arrays:
-        assert dataset_2.field_arrays[name] == \
-               pytest.approx(dataset.field_arrays[name])
+        # check data is the same
+        for attr in ('n_cells', 'n_points', 'n_arrays'):
+            if hasattr(dataset, attr):
+                assert getattr(dataset_2, attr) == getattr(dataset, attr)
+
+        for attr in ('cells', 'points'):
+            if hasattr(dataset, attr):
+                arr_have = getattr(dataset_2, attr)
+                arr_expected = getattr(dataset, attr)
+                assert arr_have == pytest.approx(arr_expected)
+
+        for name in dataset.point_arrays:
+            arr_have = dataset_2.point_arrays[name]
+            arr_expected = dataset.point_arrays[name]
+            assert arr_have == pytest.approx(arr_expected)
+
+        for name in dataset.cell_arrays:
+            arr_have = dataset_2.cell_arrays[name]
+            arr_expected = dataset.cell_arrays[name]
+            assert arr_have == pytest.approx(arr_expected)
+
+        for name in dataset.field_arrays:
+            arr_have = dataset_2.field_arrays[name]
+            arr_expected = dataset.field_arrays[name]
+            assert arr_have == pytest.approx(arr_expected)
diff --git a/tests/test_filters.py b/tests/test_filters.py
index 4a808018..d8f123d2 100644
--- a/tests/test_filters.py
+++ b/tests/test_filters.py
@@ -13,17 +13,8 @@ from pyvista import examples
 from pyvista.core.errors import VTKVersionError
 
 
-DATASETS = [
-    examples.load_uniform(),  # UniformGrid
-    examples.load_rectilinear(),  # RectilinearGrid
-    examples.load_hexbeam(),  # UnstructuredGrid
-    examples.load_airplane(),  # PolyData
-    examples.load_structured(),  # StructuredGrid
-]
 normals = ['x', 'y', '-z', (1, 1, 1), (3.3, 5.4, 0.8)]
 
-COMPOSITE = pyvista.MultiBlock(DATASETS, deep=True)
-
 skip_py2_nobind = pytest.mark.skipif(int(sys.version[0]) < 3,
                                      reason="Python 2 doesn't support binding methods")
 
@@ -31,6 +22,11 @@ skip_windows = pytest.mark.skipif(os.name == 'nt', reason="Flaky Windows tests")
 skip_mac = pytest.mark.skipif(platform.system() == 'Darwin', reason="Flaky Mac tests")
 
 
+@pytest.fixture
+def composite(datasets):
+    return pyvista.MultiBlock(datasets)
+
+
 @pytest.fixture()
 def grid():
     return pyvista.UnstructuredGrid(examples.hexbeamfile)
@@ -51,9 +47,9 @@ def test_datasetfilters_init():
 
 
 @skip_windows
-def test_clip_filter():
+def test_clip_filter(datasets):
     """This tests the clip filter on all datatypes available filters"""
-    for i, dataset in enumerate(DATASETS):
+    for i, dataset in enumerate(datasets):
         clp = dataset.clip(normal=normals[i], invert=True)
         assert clp is not None
         if isinstance(dataset, pyvista.PolyData):
@@ -62,7 +58,7 @@ def test_clip_filter():
             assert isinstance(clp, pyvista.UnstructuredGrid)
 
     # clip with get_clipped=True
-    for i, dataset in enumerate(DATASETS):
+    for i, dataset in enumerate(datasets):
         clp1, clp2 = dataset.clip(normal=normals[i], invert=True, return_clipped=True)
         for clp in (clp1, clp2):
             if isinstance(dataset, pyvista.PolyData):
@@ -72,9 +68,9 @@ def test_clip_filter():
 
 @skip_windows
 @skip_mac
-def test_clip_by_scalars_filter():
+def test_clip_by_scalars_filter(datasets):
     """This tests the clip filter on all datatypes available filters"""
-    for i, dataset_in in enumerate(DATASETS):
+    for i, dataset_in in enumerate(datasets):
         dataset = dataset_in.copy()  # don't modify in-place
         if dataset.active_scalars_info.name is None:
             dataset['scalars'] = np.arange(dataset.n_points)
@@ -91,14 +87,14 @@ def test_clip_by_scalars_filter():
 
 
 @skip_py2_nobind
-def test_clip_filter_composite():
+def test_clip_filter_composite(composite):
     # Now test composite data structures
-    output = COMPOSITE.clip(normal=normals[0], invert=False)
-    assert output.n_blocks == COMPOSITE.n_blocks
+    output = composite.clip(normal=normals[0], invert=False)
+    assert output.n_blocks == composite.n_blocks
 
 
-def test_clip_box():
-    for i, dataset in enumerate(DATASETS):
+def test_clip_box(datasets):
+    for i, dataset in enumerate(datasets):
         clp = dataset.clip_box(invert=True)
         assert clp is not None
         assert isinstance(clp, pyvista.UnstructuredGrid)
@@ -128,10 +124,10 @@ def test_clip_box():
 
 
 @skip_py2_nobind
-def test_clip_box_composite():
+def test_clip_box_composite(composite):
     # Now test composite data structures
-    output = COMPOSITE.clip_box(invert=False)
-    assert output.n_blocks == COMPOSITE.n_blocks
+    output = composite.clip_box(invert=False)
+    assert output.n_blocks == composite.n_blocks
 
 
 def test_clip_surface():
@@ -169,9 +165,9 @@ def test_implicit_distance():
     assert "implicit_distance" in dataset.point_arrays
 
 
-def test_slice_filter():
+def test_slice_filter(datasets):
     """This tests the slice filter on all datatypes available filters"""
-    for i, dataset in enumerate(DATASETS):
+    for i, dataset in enumerate(datasets):
         slc = dataset.slice(normal=normals[i])
         assert slc is not None
         assert isinstance(slc, pyvista.PolyData)
@@ -184,16 +180,16 @@ def test_slice_filter():
 
 
 @skip_py2_nobind
-def test_slice_filter_composite():
+def test_slice_filter_composite(composite):
     # Now test composite data structures
-    output = COMPOSITE.slice(normal=normals[0])
-    assert output.n_blocks == COMPOSITE.n_blocks
+    output = composite.slice(normal=normals[0])
+    assert output.n_blocks == composite.n_blocks
 
 
-def test_slice_orthogonal_filter():
+def test_slice_orthogonal_filter(datasets):
     """This tests the slice filter on all datatypes available filters"""
 
-    for i, dataset in enumerate(DATASETS):
+    for i, dataset in enumerate(datasets):
         slices = dataset.slice_orthogonal()
         assert slices is not None
         assert isinstance(slices, pyvista.MultiBlock)
@@ -203,17 +199,17 @@ def test_slice_orthogonal_filter():
 
 
 @skip_py2_nobind
-def test_slice_orthogonal_filter_composite():
+def test_slice_orthogonal_filter_composite(composite):
     # Now test composite data structures
-    output = COMPOSITE.slice_orthogonal()
-    assert output.n_blocks == COMPOSITE.n_blocks
+    output = composite.slice_orthogonal()
+    assert output.n_blocks == composite.n_blocks
 
 
-def test_slice_along_axis():
+def test_slice_along_axis(datasets):
     """Test the many slices along axis filter """
     axii = ['x', 'y', 'z', 'y', 0]
     ns = [2, 3, 4, 10, 20, 13]
-    for i, dataset in enumerate(DATASETS):
+    for i, dataset in enumerate(datasets):
         slices = dataset.slice_along_axis(n=ns[i], axis=axii[i])
         assert slices is not None
         assert isinstance(slices, pyvista.MultiBlock)
@@ -226,14 +222,14 @@ def test_slice_along_axis():
 
 
 @skip_py2_nobind
-def test_slice_along_axis_composite():
+def test_slice_along_axis_composite(composite):
     # Now test composite data structures
-    output = COMPOSITE.slice_along_axis()
-    assert output.n_blocks == COMPOSITE.n_blocks
+    output = composite.slice_along_axis()
+    assert output.n_blocks == composite.n_blocks
 
 
-def test_threshold():
-    for i, dataset in enumerate(DATASETS[0:3]):
+def test_threshold(datasets):
+    for i, dataset in enumerate(datasets[0:3]):
         thresh = dataset.threshold()
         assert thresh is not None
         assert isinstance(thresh, pyvista.UnstructuredGrid)
@@ -253,7 +249,7 @@ def test_threshold():
         dataset.threshold({100, 500})
     # Now test DATASETS without arrays
     with pytest.raises(ValueError):
-        for i, dataset in enumerate(DATASETS[3:-1]):
+        for i, dataset in enumerate(datasets[3:-1]):
             thresh = dataset.threshold()
             assert thresh is not None
             assert isinstance(thresh, pyvista.UnstructuredGrid)
@@ -261,14 +257,14 @@ def test_threshold():
     with pytest.raises(ValueError):
         dataset.threshold([10, 100, 300])
     with pytest.raises(ValueError):
-        DATASETS[0].threshold([10, 500], scalars='Spatial Point Data',
+        datasets[0].threshold([10, 500], scalars='Spatial Point Data',
                               all_scalars=True)
 
-def test_threshold_percent():
+def test_threshold_percent(datasets):
     percents = [25, 50, [18.0, 85.0], [19.0, 80.0], 0.70]
     inverts = [False, True, False, True, False]
     # Only test data sets that have arrays
-    for i, dataset in enumerate(DATASETS[0:3]):
+    for i, dataset in enumerate(datasets[0:3]):
         thresh = dataset.threshold_percent(percent=percents[i], invert=inverts[i])
         assert thresh is not None
         assert isinstance(thresh, pyvista.UnstructuredGrid)
@@ -283,68 +279,68 @@ def test_threshold_percent():
         dataset.threshold_percent({18.0, 85.0})
 
 
-def test_outline():
-    for i, dataset in enumerate(DATASETS):
+def test_outline(datasets):
+    for i, dataset in enumerate(datasets):
         outline = dataset.outline()
         assert outline is not None
         assert isinstance(outline, pyvista.PolyData)
 
 
 @skip_py2_nobind
-def test_outline_composite():
+def test_outline_composite(composite):
     # Now test composite data structures
-    output = COMPOSITE.outline()
+    output = composite.outline()
     assert isinstance(output, pyvista.PolyData)
-    output = COMPOSITE.outline(nested=True)
+    output = composite.outline(nested=True)
 
     # vtk 9.0.0 returns polydata
     assert isinstance(output, (pyvista.MultiBlock, pyvista.PolyData))
     if isinstance(output, pyvista.MultiBlock):
-        assert output.n_blocks == COMPOSITE.n_blocks
+        assert output.n_blocks == composite.n_blocks
 
 
-def test_outline_corners():
-    for i, dataset in enumerate(DATASETS):
+def test_outline_corners(datasets):
+    for i, dataset in enumerate(datasets):
         outline = dataset.outline_corners()
         assert outline is not None
         assert isinstance(outline, pyvista.PolyData)
 
 
 @skip_py2_nobind
-def test_outline_corners_composite():
+def test_outline_corners_composite(composite):
     # Now test composite data structures
-    output = COMPOSITE.outline_corners()
+    output = composite.outline_corners()
     assert isinstance(output, pyvista.PolyData)
-    output = COMPOSITE.outline_corners(nested=True)
-    assert output.n_blocks == COMPOSITE.n_blocks
+    output = composite.outline_corners(nested=True)
+    assert output.n_blocks == composite.n_blocks
 
 
-def test_extract_geometry():
-    for i, dataset in enumerate(DATASETS):
+def test_extract_geometry(datasets, composite):
+    for i, dataset in enumerate(datasets):
         outline = dataset.extract_geometry()
         assert outline is not None
         assert isinstance(outline, pyvista.PolyData)
     # Now test composite data structures
-    output = COMPOSITE.extract_geometry()
+    output = composite.extract_geometry()
     assert isinstance(output, pyvista.PolyData)
 
 
-def test_wireframe():
-    for i, dataset in enumerate(DATASETS):
+def test_wireframe(datasets):
+    for i, dataset in enumerate(datasets):
         wire = dataset.extract_all_edges()
         assert wire is not None
         assert isinstance(wire, pyvista.PolyData)
 
 
 @skip_py2_nobind
-def test_wireframe_composite():
+def test_wireframe_composite(composite):
     # Now test composite data structures
-    output = COMPOSITE.extract_all_edges()
-    assert output.n_blocks == COMPOSITE.n_blocks
+    output = composite.extract_all_edges()
+    assert output.n_blocks == composite.n_blocks
 
 
-def test_delaunay_2d():
-    mesh = DATASETS[2].delaunay_2d()  # UnstructuredGrid
+def test_delaunay_2d(datasets):
+    mesh = datasets[2].delaunay_2d()  # UnstructuredGrid
     assert isinstance(mesh, pyvista.PolyData)
     assert mesh.n_points
 
@@ -411,10 +407,10 @@ def test_elevation():
 
 
 @skip_py2_nobind
-def test_elevation_composite():
+def test_elevation_composite(composite):
     # Now test composite data structures
-    output = COMPOSITE.elevation()
-    assert output.n_blocks == COMPOSITE.n_blocks
+    output = composite.elevation()
+    assert output.n_blocks == composite.n_blocks
 
 
 def test_texture_map_to_plane():
@@ -449,8 +445,8 @@ def test_texture_map_to_sphere():
     assert 'Texture Coordinates' in dataset.array_names
 
 
-def test_compute_cell_sizes():
-    for i, dataset in enumerate(DATASETS):
+def test_compute_cell_sizes(datasets):
+    for i, dataset in enumerate(datasets):
         result = dataset.compute_cell_sizes()
         assert result is not None
         assert isinstance(result, type(dataset))
@@ -463,59 +459,55 @@ def test_compute_cell_sizes():
 
 
 @skip_py2_nobind
-def test_compute_cell_sizes_composite():
+def test_compute_cell_sizes_composite(composite):
     # Now test composite data structures
-    output = COMPOSITE.compute_cell_sizes()
-    assert output.n_blocks == COMPOSITE.n_blocks
+    output = composite.compute_cell_sizes()
+    assert output.n_blocks == composite.n_blocks
 
 
-def test_cell_centers():
-    for i, dataset in enumerate(DATASETS):
+def test_cell_centers(datasets):
+    for i, dataset in enumerate(datasets):
         result = dataset.cell_centers()
         assert result is not None
         assert isinstance(result, pyvista.PolyData)
 
 
 @skip_py2_nobind
-def test_cell_centers_composite():
+def test_cell_centers_composite(composite):
     # Now test composite data structures
-    output = COMPOSITE.cell_centers()
-    assert output.n_blocks == COMPOSITE.n_blocks
+    output = composite.cell_centers()
+    assert output.n_blocks == composite.n_blocks
 
 
-def test_glyph():
-    for i, dataset in enumerate(DATASETS):
+def test_glyph(datasets, sphere):
+    for i, dataset in enumerate(datasets):
+        dataset.vectors = np.ones_like(dataset.points)
         result = dataset.glyph()
         assert result is not None
         assert isinstance(result, pyvista.PolyData)
     # Test different options for glyph filter
-    sphere = pyvista.Sphere(radius=3.14, theta_resolution=5, phi_resolution=5)
     sphere_sans_arrays = sphere.copy()
-    # make cool swirly pattern
-    vectors = np.vstack((np.sin(sphere.points[:, 0]),
-                        np.cos(sphere.points[:, 1]),
-                        np.cos(sphere.points[:, 2]))).T
-    # add and scale
-    sphere.vectors = vectors*0.3
-    sphere.point_arrays['foo'] = np.random.rand(sphere.n_points)
+    sphere.compute_normals(inplace=True)
+    sphere.vectors = np.ones([sphere.n_points,3])
     sphere.point_arrays['arr'] = np.ones(sphere.n_points)
-    result = sphere.glyph(scale=False)
-    result = sphere.glyph(scale='arr')
-    result = sphere.glyph(scale='arr', orient='Normals', factor=0.1)
-    result = sphere.glyph(scale='arr', orient='Normals', factor=0.1, tolerance=0.1)
-    result = sphere.glyph(scale='arr', orient='Normals', factor=0.1, tolerance=0.1,
+
+    assert sphere.glyph(scale=False)
+    assert sphere.glyph(scale='arr')
+    assert sphere.glyph(scale='arr', orient='Normals', factor=0.1)
+    assert sphere.glyph(scale='arr', orient='Normals', factor=0.1, tolerance=0.1)
+    assert sphere.glyph(scale='arr', orient='Normals', factor=0.1, tolerance=0.1,
                           clamping=False, rng=[1, 1])
     # passing one or more custom glyphs; many cases for full coverage
     geoms = [pyvista.Sphere(theta_resolution=5, phi_resolution=5),
              pyvista.Arrow(tip_resolution=5, shaft_resolution=5),
              pyvista.ParametricSuperToroid(u_res=10, v_res=10, w_res=10)]
     indices = range(len(geoms))
-    result = sphere.glyph(geom=geoms[0])
-    result = sphere.glyph(geom=geoms, indices=indices, rng=(0, len(geoms)))
-    result = sphere.glyph(geom=geoms)
-    result = sphere.glyph(geom=geoms, scale='arr', orient='Normals', factor=0.1, tolerance=0.1)
-    result = sphere.glyph(geom=geoms[:1], indices=[None])
-    result = sphere_sans_arrays.glyph(geom=geoms)
+    assert sphere.glyph(geom=geoms[0])
+    assert sphere.glyph(geom=geoms, indices=indices, rng=(0, len(geoms)))
+    assert sphere.glyph(geom=geoms)
+    assert sphere.glyph(geom=geoms, scale='arr', orient='Normals', factor=0.1, tolerance=0.1)
+    assert sphere.glyph(geom=geoms[:1], indices=[None])
+    assert sphere_sans_arrays.glyph(geom=geoms)
     with pytest.raises(TypeError):
         # wrong type for the glyph
         sphere.glyph(geom=pyvista.StructuredGrid())
@@ -527,6 +519,20 @@ def test_glyph():
         sphere.glyph(geom=geoms, indices=indices[:-1])
 
 
+def test_glyph_cell_point_data(sphere):
+    sphere['vectors_cell'] = np.ones([sphere.n_cells,3])
+    sphere['vectors_points'] = np.ones([sphere.n_points,3])
+    sphere['arr_cell'] = np.ones(sphere.n_cells)
+    sphere['arr_points'] = np.ones(sphere.n_points)
+    
+    assert sphere.glyph(orient='vectors_cell', scale='arr_cell')
+    assert sphere.glyph(orient='vectors_points', scale='arr_points')
+    with pytest.raises(ValueError):
+        sphere.glyph(orient='vectors_cell', scale='arr_points')
+    with pytest.raises(ValueError):
+        sphere.glyph(orient='vectors_points', scale='arr_cell')
+
+
 def test_split_and_connectivity():
     # Load a simple example mesh
     dataset = examples.load_uniform()
@@ -605,10 +611,10 @@ def test_cell_data_to_point_data():
 
 
 @skip_py2_nobind
-def test_cell_data_to_point_data_composite():
+def test_cell_data_to_point_data_composite(composite):
     # Now test composite data structures
-    output = COMPOSITE.cell_data_to_point_data()
-    assert output.n_blocks == COMPOSITE.n_blocks
+    output = composite.cell_data_to_point_data()
+    assert output.n_blocks == composite.n_blocks
 
 
 def test_point_data_to_cell_data():
@@ -620,10 +626,10 @@ def test_point_data_to_cell_data():
 
 
 @skip_py2_nobind
-def test_point_data_to_cell_data_composite():
+def test_point_data_to_cell_data_composite(composite):
     # Now test composite data structures
-    output = COMPOSITE.point_data_to_cell_data()
-    assert output.n_blocks == COMPOSITE.n_blocks
+    output = composite.point_data_to_cell_data()
+    assert output.n_blocks == composite.n_blocks
 
 
 def test_triangulate():
@@ -634,10 +640,10 @@ def test_triangulate():
 
 
 @skip_py2_nobind
-def test_triangulate_composite():
+def test_triangulate_composite(composite):
     # Now test composite data structures
-    output = COMPOSITE.triangulate()
-    assert output.n_blocks == COMPOSITE.n_blocks
+    output = composite.triangulate()
+    assert output.n_blocks == composite.n_blocks
 
 
 def test_delaunay_3d():
@@ -719,6 +725,7 @@ def test_streamlines_start_position(uniform_vec):
 
     assert all([stream.n_points, stream.n_cells])
 
+
 def test_streamlines_errors(uniform_vec):
     with pytest.raises(ValueError):
         uniform_vec.streamlines('vectors', integration_direction='not valid')
@@ -983,14 +990,15 @@ def test_extract_points():
     assert sub_surf_adj.n_cells == 4
     assert sub_surf_nocells.cells[0] == 1
 
+
 @skip_py2_nobind
-def test_slice_along_line_composite():
+def test_slice_along_line_composite(composite):
     # Now test composite data structures
-    a = [COMPOSITE.bounds[0], COMPOSITE.bounds[2], COMPOSITE.bounds[4]]
-    b = [COMPOSITE.bounds[1], COMPOSITE.bounds[3], COMPOSITE.bounds[5]]
+    a = [composite.bounds[0], composite.bounds[2], composite.bounds[4]]
+    b = [composite.bounds[1], composite.bounds[3], composite.bounds[5]]
     line = pyvista.Line(a, b, resolution=10)
-    output = COMPOSITE.slice_along_line(line)
-    assert output.n_blocks == COMPOSITE.n_blocks
+    output = composite.slice_along_line(line)
+    assert output.n_blocks == composite.n_blocks
 
 
 def test_interpolate():
@@ -1287,71 +1295,73 @@ def test_shrink():
     assert shrunk.area < mesh.area
 
 
-@pytest.mark.parametrize('dataset,num_cell_arrays,num_point_arrays',
-                         itertools.product(DATASETS, [0, 1, 2], [0, 1, 2]))
-def test_transform_mesh(dataset, num_cell_arrays, num_point_arrays):
+@pytest.mark.parametrize('num_cell_arrays,num_point_arrays',
+                         itertools.product([0, 1, 2], [0, 1, 2]))
+def test_transform_mesh(datasets, num_cell_arrays, num_point_arrays):
     # rotate about x-axis by 90 degrees
-    tf = pyvista.transformations.axis_angle_rotation((1, 0, 0), 90)
+    for dataset in datasets:
+        tf = pyvista.transformations.axis_angle_rotation((1, 0, 0), 90)
 
-    for i in range(num_cell_arrays):
-        dataset.cell_arrays['C%d' % i] = np.random.rand(dataset.n_cells, 3)
+        for i in range(num_cell_arrays):
+            dataset.cell_arrays['C%d' % i] = np.random.rand(dataset.n_cells, 3)
 
-    for i in range(num_point_arrays):
-        dataset.point_arrays['P%d' % i] = np.random.rand(dataset.n_points, 3)
+        for i in range(num_point_arrays):
+            dataset.point_arrays['P%d' % i] = np.random.rand(dataset.n_points, 3)
 
-    # deactivate any active vectors!
-    # even if transform_all_input_vectors is False, vtkTransformfilter will
-    # transform active vectors
-    dataset.set_active_vectors(None)
+        # deactivate any active vectors!
+        # even if transform_all_input_vectors is False, vtkTransformfilter will
+        # transform active vectors
+        dataset.set_active_vectors(None)
 
-    transformed = dataset.transform(tf, transform_all_input_vectors=False, inplace=False)
+        transformed = dataset.transform(tf, transform_all_input_vectors=False, inplace=False)
 
-    assert dataset.points[:, 0] == pytest.approx(transformed.points[:, 0])
-    assert dataset.points[:, 2] == pytest.approx(-transformed.points[:, 1])
-    assert dataset.points[:, 1] == pytest.approx(transformed.points[:, 2])
+        assert dataset.points[:, 0] == pytest.approx(transformed.points[:, 0])
+        assert dataset.points[:, 2] == pytest.approx(-transformed.points[:, 1])
+        assert dataset.points[:, 1] == pytest.approx(transformed.points[:, 2])
 
-    # ensure that none of the vector data is changed
-    for name, array in dataset.point_arrays.items():
-        assert transformed.point_arrays[name] == pytest.approx(array)
+        # ensure that none of the vector data is changed
+        for name, array in dataset.point_arrays.items():
+            assert transformed.point_arrays[name] == pytest.approx(array)
 
-    for name, array in dataset.cell_arrays.items():
-        assert transformed.cell_arrays[name] == pytest.approx(array)
+        for name, array in dataset.cell_arrays.items():
+            assert transformed.cell_arrays[name] == pytest.approx(array)
 
 
-@pytest.mark.parametrize('dataset,num_cell_arrays,num_point_arrays',
-                         itertools.product(DATASETS, [0, 1, 2], [0, 1, 2]))
-def test_transform_mesh_and_vectors(dataset, num_cell_arrays, num_point_arrays):
-      # rotate about x-axis by 90 degrees
-    tf = pyvista.transformations.axis_angle_rotation((1, 0, 0), 90)
+@pytest.mark.parametrize('num_cell_arrays,num_point_arrays',
+                         itertools.product([0, 1, 2], [0, 1, 2]))
+def test_transform_mesh_and_vectors(datasets, num_cell_arrays, num_point_arrays):
+    for dataset in datasets:
+        # rotate about x-axis by 90 degrees
+        tf = pyvista.transformations.axis_angle_rotation((1, 0, 0), 90)
 
-    for i in range(num_cell_arrays):
-        dataset.cell_arrays['C%d' % i] = np.random.rand(dataset.n_cells, 3)
+        for i in range(num_cell_arrays):
+            dataset.cell_arrays['C%d' % i] = np.random.rand(dataset.n_cells, 3)
 
-    for i in range(num_point_arrays):
-        dataset.point_arrays['P%d' % i] = np.random.rand(dataset.n_points, 3)
+        for i in range(num_point_arrays):
+            dataset.point_arrays['P%d' % i] = np.random.rand(dataset.n_points, 3)
 
-    # handle
-    f = pyvista._vtk.vtkTransformFilter()
-    if not hasattr(f, 'SetTransformAllInputVectors'):
-        with pytest.raises(VTKVersionError):
-            transformed = dataset.transform(tf, transform_all_input_vectors=True, inplace=False)
-        return
+        # handle
+        f = pyvista._vtk.vtkTransformFilter()
+        if not hasattr(f, 'SetTransformAllInputVectors'):
+            with pytest.raises(VTKVersionError):
+                transformed = dataset.transform(tf, transform_all_input_vectors=True, inplace=False)
+            return
 
-    transformed = dataset.transform(tf, transform_all_input_vectors=True, inplace=False)
+        transformed = dataset.transform(tf, transform_all_input_vectors=True, inplace=False)
 
-    assert dataset.points[:, 0] == pytest.approx(transformed.points[:, 0])
-    assert dataset.points[:, 2] == pytest.approx(-transformed.points[:, 1])
-    assert dataset.points[:, 1] == pytest.approx(transformed.points[:, 2])
+        assert dataset.points[:, 0] == pytest.approx(transformed.points[:, 0])
+        assert dataset.points[:, 2] == pytest.approx(-transformed.points[:, 1])
+        assert dataset.points[:, 1] == pytest.approx(transformed.points[:, 2])
 
-    for i in range(num_cell_arrays):
-        assert dataset.cell_arrays['C%d' % i][:, 0] == pytest.approx( transformed.cell_arrays['C%d' % i][:, 0])
-        assert dataset.cell_arrays['C%d' % i][:, 2] == pytest.approx(-transformed.cell_arrays['C%d' % i][:, 1])
-        assert dataset.cell_arrays['C%d' % i][:, 1] == pytest.approx( transformed.cell_arrays['C%d' % i][:, 2])
+        for i in range(num_cell_arrays):
+            assert dataset.cell_arrays['C%d' % i][:, 0] == pytest.approx( transformed.cell_arrays['C%d' % i][:, 0])
+            assert dataset.cell_arrays['C%d' % i][:, 2] == pytest.approx(-transformed.cell_arrays['C%d' % i][:, 1])
+            assert dataset.cell_arrays['C%d' % i][:, 1] == pytest.approx( transformed.cell_arrays['C%d' % i][:, 2])
 
-    for i in range(num_point_arrays):
-        assert dataset.point_arrays['P%d' % i][:, 0] == pytest.approx( transformed.point_arrays['P%d' % i][:, 0])
-        assert dataset.point_arrays['P%d' % i][:, 2] == pytest.approx(-transformed.point_arrays['P%d' % i][:, 1])
-        assert dataset.point_arrays['P%d' % i][:, 1] == pytest.approx( transformed.point_arrays['P%d' % i][:, 2])
+        for i in range(num_point_arrays):
+            assert dataset.point_arrays['P%d' % i][:, 0] == pytest.approx( transformed.point_arrays['P%d' % i][:, 0])
+            assert dataset.point_arrays['P%d' % i][:, 2] == pytest.approx(-transformed.point_arrays['P%d' % i][:, 1])
+            assert dataset.point_arrays['P%d' % i][:, 1] == pytest.approx( transformed.point_arrays['P%d' % i][:, 2])
 
 
 @pytest.mark.parametrize('dataset', [
@@ -1365,48 +1375,48 @@ def test_transform_inplace_bad_types(dataset):
         dataset.transform(tf, inplace=True)
 
 
-@pytest.mark.parametrize('dataset', DATASETS)
-def test_reflect_mesh_about_point(dataset):
-    x_plane = 500
-    reflected = dataset.reflect((1, 0, 0), point=(x_plane, 0, 0))
-    assert reflected.n_cells == dataset.n_cells
-    assert reflected.n_points == dataset.n_points
-    assert np.allclose(x_plane - dataset.points[:, 0], reflected.points[:, 0] - x_plane)
-    assert np.allclose(dataset.points[:, 1:], reflected.points[:, 1:])
+def test_reflect_mesh_about_point(datasets):
+    for dataset in datasets:
+        x_plane = 500
+        reflected = dataset.reflect((1, 0, 0), point=(x_plane, 0, 0))
+        assert reflected.n_cells == dataset.n_cells
+        assert reflected.n_points == dataset.n_points
+        assert np.allclose(x_plane - dataset.points[:, 0], reflected.points[:, 0] - x_plane)
+        assert np.allclose(dataset.points[:, 1:], reflected.points[:, 1:])
 
 
 @pytest.mark.skipif(not VTK9, reason='Only supported on VTK v9 or newer')
-@pytest.mark.parametrize('dataset', DATASETS)
-def test_reflect_mesh_with_vectors(dataset):
-    if hasattr(dataset, 'compute_normals'):
-        dataset.compute_normals(inplace=True)
-
-    # add vector data to cell and point arrays
-    dataset.cell_arrays['C'] = np.arange(dataset.n_cells)[:, np.newaxis] * \
-                                np.array([1, 2, 3], dtype=float).reshape((1, 3))
-    dataset.point_arrays['P'] = np.arange(dataset.n_points)[:, np.newaxis] * \
-                                np.array([1, 2, 3], dtype=float).reshape((1, 3))
-
-    reflected = dataset.reflect((1, 0, 0), transform_all_input_vectors=True, inplace=False)
-
-    # assert isinstance(reflected, type(dataset))
-    assert reflected.n_cells == dataset.n_cells
-    assert reflected.n_points == dataset.n_points
-    assert np.allclose(dataset.points[:, 0], -reflected.points[:, 0])
-    assert np.allclose(dataset.points[:, 1:], reflected.points[:, 1:])
-
-    # assert normals are reflected
-    if hasattr(dataset, 'compute_normals'):
-        assert np.allclose(dataset.cell_arrays['Normals'][:, 0], -reflected.cell_arrays['Normals'][:, 0])
-        assert np.allclose(dataset.cell_arrays['Normals'][:, 1:], reflected.cell_arrays['Normals'][:, 1:])
-        assert np.allclose(dataset.point_arrays['Normals'][:, 0], -reflected.point_arrays['Normals'][:, 0])
-        assert np.allclose(dataset.point_arrays['Normals'][:, 1:], reflected.point_arrays['Normals'][:, 1:])
-
-    # assert other vector fields are reflected
-    assert np.allclose(dataset.cell_arrays['C'][:, 0], -reflected.cell_arrays['C'][:, 0])
-    assert np.allclose(dataset.cell_arrays['C'][:, 1:], reflected.cell_arrays['C'][:, 1:])
-    assert np.allclose(dataset.point_arrays['P'][:, 0], -reflected.point_arrays['P'][:, 0])
-    assert np.allclose(dataset.point_arrays['P'][:, 1:], reflected.point_arrays['P'][:, 1:])
+def test_reflect_mesh_with_vectors(datasets):
+    for dataset in datasets:
+        if hasattr(dataset, 'compute_normals'):
+            dataset.compute_normals(inplace=True)
+
+        # add vector data to cell and point arrays
+        dataset.cell_arrays['C'] = np.arange(dataset.n_cells)[:, np.newaxis] * \
+            np.array([1, 2, 3], dtype=float).reshape((1, 3))
+        dataset.point_arrays['P'] = np.arange(dataset.n_points)[:, np.newaxis] * \
+            np.array([1, 2, 3], dtype=float).reshape((1, 3))
+
+        reflected = dataset.reflect((1, 0, 0), transform_all_input_vectors=True, inplace=False)
+
+        # assert isinstance(reflected, type(dataset))
+        assert reflected.n_cells == dataset.n_cells
+        assert reflected.n_points == dataset.n_points
+        assert np.allclose(dataset.points[:, 0], -reflected.points[:, 0])
+        assert np.allclose(dataset.points[:, 1:], reflected.points[:, 1:])
+
+        # assert normals are reflected
+        if hasattr(dataset, 'compute_normals'):
+            assert np.allclose(dataset.cell_arrays['Normals'][:, 0], -reflected.cell_arrays['Normals'][:, 0])
+            assert np.allclose(dataset.cell_arrays['Normals'][:, 1:], reflected.cell_arrays['Normals'][:, 1:])
+            assert np.allclose(dataset.point_arrays['Normals'][:, 0], -reflected.point_arrays['Normals'][:, 0])
+            assert np.allclose(dataset.point_arrays['Normals'][:, 1:], reflected.point_arrays['Normals'][:, 1:])
+
+        # assert other vector fields are reflected
+        assert np.allclose(dataset.cell_arrays['C'][:, 0], -reflected.cell_arrays['C'][:, 0])
+        assert np.allclose(dataset.cell_arrays['C'][:, 1:], reflected.cell_arrays['C'][:, 1:])
+        assert np.allclose(dataset.point_arrays['P'][:, 0], -reflected.point_arrays['P'][:, 0])
+        assert np.allclose(dataset.point_arrays['P'][:, 1:], reflected.point_arrays['P'][:, 1:])
 
 
 @pytest.mark.parametrize('dataset', [
diff --git a/tests/test_polydata.py b/tests/test_polydata.py
index 1f019b66..e407ca58 100644
--- a/tests/test_polydata.py
+++ b/tests/test_polydata.py
@@ -36,11 +36,6 @@ def cube_dense():
 
 test_path = os.path.dirname(os.path.abspath(__file__))
 
-try:
-    CONDA_ENV = os.environ['CONDA_ALWAYS_YES'] == "1"
-except KeyError:
-    CONDA_ENV = False
-
 
 def is_binary(filename):
     """Return ``True`` when a file is binary."""
@@ -248,7 +243,6 @@ def test_ray_trace_plot(sphere):
     assert np.any(ind)
 
 
-@pytest.mark.skipif(not CONDA_ENV, reason="Requires libspatialindex dependency only installable via conda")
 def test_multi_ray_trace(sphere):
     pytest.importorskip('rtree')
     pytest.importorskip('pyembree')