diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml
new file mode 100644
index 0000000000000000000000000000000000000000..a7d4c6d4cea39193fb02d5b06c3116f58bfd1a19
--- /dev/null
+++ b/.test-conda-env-py3.yml
@@ -0,0 +1,39 @@
+name: test-conda-env
+channels:
+- conda-forge
+- defaults
+
+dependencies:
+- python=3
+- git
+- conda-forge::numpy
+- pocl
+- mako
+- pyopencl
+- islpy
+- gmsh
+
+# for Pytential
+- cython
+
+# for pymetis
+- pybind11
+
+# Only needed to make pylint succeed
+- matplotlib
+
+- pip
+- pip:
+    - git+https://github.com/inducer/pymbolic
+    - git+https://github.com/inducer/loopy
+
+    - git+https://gitlab.tiker.net/inducer/gmsh_interop.git
+    - git+https://gitlab.tiker.net/inducer/modepy.git
+
+    # required by pytential, which is in turn needed for some tests
+    - git+https://gitlab.tiker.net/inducer/boxtree.git
+    - git+https://gitlab.tiker.net/inducer/sumpy.git
+    - git+https://gitlab.tiker.net/inducer/pytential.git
+
+    # requires pymetis for tests for partition_mesh
+    - git+https://gitlab.tiker.net/inducer/pymetis.git
diff --git a/README.rst b/README.rst
index 1a033e2e0ad3c1acbd0cef5b12cf7905fa24782d..b012a2185b14361a68562f5fb50cfbf5b6041fbb 100644
--- a/README.rst
+++ b/README.rst
@@ -2,9 +2,14 @@ meshmode: High-Order Meshes and Discontinuous Function Spaces
 =============================================================
 
 .. image:: https://gitlab.tiker.net/inducer/meshmode/badges/master/pipeline.svg
-   :target: https://gitlab.tiker.net/inducer/meshmode/commits/master
+    :alt: Gitlab Build Status
+    :target: https://gitlab.tiker.net/inducer/meshmode/commits/master
+.. image:: https://dev.azure.com/ak-spam/inducer/_apis/build/status/inducer.meshmode?branchName=master
+    :alt: Azure Build Status
+    :target: https://dev.azure.com/ak-spam/inducer/_build/latest?definitionId=14&branchName=master
 .. image:: https://badge.fury.io/py/meshmode.png
-  :target: http://pypi.python.org/pypi/meshmode
+    :alt: Python Package Index Release Page
+    :target: https://pypi.org/project/meshmode/
 
 * `Source code on Github <https://github.com/inducer/meshmode>`_
 * `Documentation <https://documen.tician.de/meshmode>`_
diff --git a/azure-pipelines.yml b/azure-pipelines.yml
new file mode 100644
index 0000000000000000000000000000000000000000..4f83fdd1874947da829e5db75d1b49571fa8f44a
--- /dev/null
+++ b/azure-pipelines.yml
@@ -0,0 +1,66 @@
+jobs:
+-
+    job: 'Python2'
+    pool:
+        vmImage: 'ubuntu-latest'
+
+    steps:
+    -
+        script: |
+            set -e
+            sed 's/python=3/python=2.7/' .test-conda-env-py3.yml > .test-conda-env-py2.yml
+            cat .test-conda-env-py2.yml
+            CONDA_ENVIRONMENT=.test-conda-env-py2.yml
+            curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project-within-miniconda.sh
+            . ./build-and-test-py-project-within-miniconda.sh
+
+        displayName: 'Pytest Conda'
+    -
+        task: PublishTestResults@2
+        inputs:
+            testResultsFormat: 'JUnit'
+            testResultsFiles: 'test/pytest.xml'
+
+-
+    job: 'Python3'
+    pool:
+        vmImage: 'ubuntu-latest'
+
+    steps:
+    -
+        script: |
+            set -e
+            CONDA_ENVIRONMENT=.test-conda-env-py3.yml
+            curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project-within-miniconda.sh
+            . ./build-and-test-py-project-within-miniconda.sh
+
+        displayName: 'Pytest Conda'
+
+    -
+        task: PublishTestResults@2
+        inputs:
+            testResultsFormat: 'JUnit'
+            testResultsFiles: 'test/pytest.xml'
+
+-
+    job: 'Flake8'
+    pool:
+        vmImage: 'ubuntu-latest'
+    strategy:
+        matrix:
+            Python37:
+                python.version: '3.7'
+
+    steps:
+    -
+        task: UsePythonVersion@0
+        inputs:
+            versionSpec: '$(python.version)'
+
+    -
+        script: |
+            set -e
+            curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/prepare-and-run-flake8.sh
+            . ./prepare-and-run-flake8.sh meshmode test
+
+        displayName: 'Flake8'
diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py
index 003bb834bdd1b206173b75f5264612687d055dc6..e985f3ffcd9f077a2c6dcb97b6e94e4c8147c49c 100644
--- a/meshmode/discretization/visualization.py
+++ b/meshmode/discretization/visualization.py
@@ -101,6 +101,7 @@ class _VisConnectivityGroup(Record):
 class Visualizer(object):
     """
     .. automethod:: show_scalar_in_mayavi
+    .. automethod:: show_scalar_in_matplotlib_3d
     .. automethod:: write_vtk_file
     """
 
@@ -348,6 +349,76 @@ class Visualizer(object):
 
         # }}}
 
+    # {{{ matplotlib 3D
+
+    def show_scalar_in_matplotlib_3d(self, field, **kwargs):
+        import matplotlib.pyplot as plt
+
+        # This import also registers the 3D projection.
+        import mpl_toolkits.mplot3d.art3d as art3d
+
+        do_show = kwargs.pop("do_show", True)
+        vmin = kwargs.pop("vmin", None)
+        vmax = kwargs.pop("vmax", None)
+        norm = kwargs.pop("norm", None)
+
+        with cl.CommandQueue(self.vis_discr.cl_context) as queue:
+            nodes = self.vis_discr.nodes().with_queue(queue).get()
+
+            field = self._resample_and_get(queue, field)
+
+        assert nodes.shape[0] == self.vis_discr.ambient_dim
+
+        vis_connectivity, = self._vis_connectivity()
+
+        fig = plt.gcf()
+        ax = fig.gca(projection="3d")
+
+        had_data = ax.has_data()
+
+        if self.vis_discr.dim == 2:
+            nodes = list(nodes)
+            # pad to 3D with zeros
+            while len(nodes) < 3:
+                nodes.append(0*nodes[0])
+
+            from matplotlib.tri.triangulation import Triangulation
+            tri, args, kwargs = \
+                Triangulation.get_from_args_and_kwargs(
+                        *nodes,
+                        triangles=vis_connectivity.vis_connectivity.reshape(-1, 3))
+
+            triangles = tri.get_masked_triangles()
+            xt = nodes[0][triangles]
+            yt = nodes[1][triangles]
+            zt = nodes[2][triangles]
+            verts = np.stack((xt, yt, zt), axis=-1)
+
+            fieldt = field[triangles]
+
+            polyc = art3d.Poly3DCollection(verts, **kwargs)
+
+            # average over the three points of each triangle
+            avg_field = fieldt.mean(axis=1)
+            polyc.set_array(avg_field)
+
+            if vmin is not None or vmax is not None:
+                polyc.set_clim(vmin, vmax)
+            if norm is not None:
+                polyc.set_norm(norm)
+
+            ax.add_collection(polyc)
+            ax.auto_scale_xyz(xt, yt, zt, had_data)
+
+        else:
+            raise RuntimeError("meshes of bulk dimension %d are currently "
+                    "unsupported" % self.vis_discr.dim)
+
+        if do_show:
+            plt.show()
+
+    # }}}
+
 
 def make_visualizer(queue, discr, vis_order, element_shrink_factor=None):
     from meshmode.discretization import Discretization
diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index ddecc69fda5ff2389ee5bbf14e5f975f31264c82..80289cacd4f21586e1e048c5a221d76195bd6a87 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -427,17 +427,26 @@ def generate_icosahedron(r, order):
 
 # {{{ generate_icosphere
 
-def generate_icosphere(r, order):
+def generate_icosphere(r, order, uniform_refinement_rounds=0):
     mesh = generate_icosahedron(r, order)
 
-    grp, = mesh.groups
+    if uniform_refinement_rounds:
+        # These come out conformal, so we're OK to use the faster refiner.
+        from meshmode.mesh.refinement import RefinerWithoutAdjacency
+        refiner = RefinerWithoutAdjacency(mesh)
+        for i in range(uniform_refinement_rounds):
+            refiner.refine_uniformly()
+
+        mesh = refiner.get_current_mesh()
 
+    vertices = mesh.vertices * r / np.sqrt(np.sum(mesh.vertices**2, axis=0))
+    grp, = mesh.groups
     grp = grp.copy(
             nodes=grp.nodes * r / np.sqrt(np.sum(grp.nodes**2, axis=0)))
 
     from meshmode.mesh import Mesh
     return Mesh(
-            mesh.vertices, [grp],
+            vertices, [grp],
             is_conforming=True)
 
 # }}}