From 897e6dfde96960257e9a29aa0ec9aca386ad9b01 Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Fri, 25 Oct 2024 16:12:04 +0300
Subject: [PATCH 1/6] start porting extended documentation to Sphinx

---
 docs/advanced.rst |   2 +
 docs/extended.rst | 635 ++++++++++++++++++++++++++++++++++++++++++++++
 docs/index.rst    |   1 +
 3 files changed, 638 insertions(+)
 create mode 100644 docs/extended.rst

diff --git a/docs/advanced.rst b/docs/advanced.rst
index d6986de3..7b59a39d 100644
--- a/docs/advanced.rst
+++ b/docs/advanced.rst
@@ -242,3 +242,5 @@ The remaining DOFs are internal to the element and not shared:
    
 Each DOF is associated either with a node (``nodal_dofs``), a facet
 (``facet_dofs``), an edge (``edge_dofs``), or an element (``interior_dofs``).
+
+
diff --git a/docs/extended.rst b/docs/extended.rst
new file mode 100644
index 00000000..bd106489
--- /dev/null
+++ b/docs/extended.rst
@@ -0,0 +1,635 @@
+.. _extended:
+
+========================
+ Extended documentation
+========================
+
+It is extremely helpful when working with ``skfem`` to understand the
+concepts of quadrature points and degrees-of-freedom, and how these
+concepts relate to symbolic math expressions, Python functions, and
+real-world data. This all happens in the umbrella concept of a "mesh"
+and which also has to be reasonably well understood to use ``skfem``. To
+illustrate how these components fit together, we start with a simple
+2D mesh of triangles created directly from ``skfem``.
+
+.. doctest::
+
+   >>> import matplotlib.pyplot as plt
+   >>> import numpy as np
+   >>> import skfem
+   >>> mesh = skfem.MeshTri()
+
+The mesh is represented as two lists of information. The first is a
+list of points describing where vertexes are and the second is a list
+of connections between the points to form elements. These connections
+are the "topology" of the mesh.
+
+.. doctest::
+
+   >>> mesh.p  # the list of points
+   array([[0., 1., 0., 1.],
+          [0., 0., 1., 1.]])
+
+
+.. doctest::
+
+   >>> mesh.t  # the topology
+   array([[0, 1],
+          [1, 2],
+          [2, 3]], dtype=int32)
+
+This is easy enough to draw in ``matplotlib``. The first row in the points
+array contains the x locations of each point and the second row
+contains the y locations.
+
+.. doctest::
+
+   >>> plt.plot(mesh.p[0], mesh.p[1], 'ok')
+
+.. plot::
+
+   import skfem as fem
+   mesh = fem.MeshTri()
+   plt.plot(mesh.p[0], mesh.p[1], 'ok')
+
+The topology (``mesh.t``) specifies how these points are connected
+together to form triangles. Each column specifes the indexes of the
+points that defined the vertexes of the triangle. We can add these
+lines to the plot.
+
+.. doctest::
+
+   >>> plt.plot(mesh.p[0], mesh.p[1], 'ok')
+   >>> for t in mesh.t.T: # transpose to iterate over columns
+   >>>     plt.plot(mesh.p[0,[t[0],t[1]]], mesh.p[1,[t[0],t[1]]], 'k')  # from vertex 0 to 1
+   >>>     plt.plot(mesh.p[0,[t[1],t[2]]], mesh.p[1,[t[1],t[2]]], 'k')  # from vertex 1 to 2
+   >>>     plt.plot(mesh.p[0,[t[2],t[0]]], mesh.p[1,[t[2],t[0]]], 'k')  # from vertex 2 back to 0
+
+.. plot::
+
+   import skfem as fem
+   mesh = fem.MeshTri()
+   plt.plot(mesh.p[0], mesh.p[1], 'ok')
+   for t in mesh.t.T: # transpose to iterate over columns
+       plt.plot(mesh.p[0,[t[0],t[1]]], mesh.p[1,[t[0],t[1]]], 'k')  # from vertex 0 to 1
+       plt.plot(mesh.p[0,[t[1],t[2]]], mesh.p[1,[t[1],t[2]]], 'k')  # from vertex 1 to 2
+       plt.plot(mesh.p[0,[t[2],t[0]]], mesh.p[1,[t[2],t[0]]], 'k')  # from vertex 2 back to 0
+
+This is the most basic triangulation of the unit square. ``skfem`` has
+utilities to refined this mesh into more and smaller triangles as well
+as translate it or map it into different coordinates. It also has
+basic visualization functions to make drawing the mesh on an ``matplotlib`` axis
+easier.
+
+.. doctest::
+
+   >>> import skfem.visuals.matplotlib
+   >>> mesh = skfem.MeshTri().refined(1)
+   >>> plt.subplots(figsize=(5,5))
+   >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())  # gca: "get current axis"
+
+.. plot::
+
+   import skfem
+   import skfem.visuals.matplotlib
+   mesh = skfem.MeshTri().refined(1)
+   plt.subplots(figsize=(5,5))
+   skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())  # gca: "get current axis"
+
+The ``skfem`` documentation and code uses several terms when working
+with meshes of one, two, or three dimensions that are worth clarifying
+before we proceed. These are cells/elements, ``facets``,
+``edges``, and ``nodes``, and are best illustrated with a picture:
+
+.. figure:: https://user-images.githubusercontent.com/38136423/144346451-e43fa714-2e12-4b31-a809-38359c9110aa.png
+
+   The naming conventions used in ``skfem``.
+
+Using this naming convention, ``facets`` are always shared between
+cells/elements and one dimension lower than the mesh. ``nodes`` are always at
+the vertices of the mesh. This picture also illustrates quadrilateral
+meshes, which are an alternative to triangulations that can be
+generated by ``skfem``. For the remainder of this discussion, we will work
+with 2D triangular meshes.
+
+Meshes form a kind of coordinate system to work in, and we construct a
+set of basis functions in this system by specifying a functional form
+over one cell/element of the mesh. This discussion will be limited to two
+kinds of basis functions: ones that are constant over the cell/element and
+ones that are linear over the cell/element. ``skfem`` calls these ``ElementTriP0`` and
+``ElementTriP1``, respectively. Note that these two basis sets have
+different continuity characteristics between cells/elements. Basis functions in
+``ElementTriP0`` are discontinuous between cells/elements. Basis functions in
+``ElementTriP1`` are continuous between adjacent cells/elements, but their
+derivatives are not.
+
+We continue this discussion by building a set of basis functions using
+``ElementTriP1`` over the once refined triangulation of the unit square
+discussed above.
+
+.. doctest::
+
+   >>> basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
+   >>> print(type(basis_p1))
+   <class 'skfem.assembly.basis.cell_basis.CellBasis'>
+
+What we get back after this call is a Python object of type
+``CellBasis``. This is a mostly opaque object that we can use to work
+with the set of basis functions that span our finite element
+space. Functions represented in this finite space are (obviously)
+described by a finite number of parameters, in ``skfem`` called the
+degrees-of-freedom (dofs). In our P1 space that we've constructed,
+this will always be equal to the number of nodes in the mesh. However,
+this is in general not true, so to get a ``numpy`` array of the correct
+length and initialized to zeros, we will use our basis object.
+
+.. doctest::
+
+   >>> fe_approximation = basis_p1.zeros()
+   
+Although this is a simple ``numpy`` array, there are not many things we
+can do with it directly, since out at this level of the code we don't
+know anything about what the array index means. Its primary
+application in our code will be controlling Dirichlet boundary
+conditions: those locations on the mesh where we already know the
+value of the solution. We can experiment with this by projecting a
+constant function of 1 into the finite element space, and then showing
+how we can manipulate this function using our ``basis_p1`` object and the
+``fe_approximation`` array. For now, we will also make use of another
+helper function from ``skfem`` to visualize the functions we
+construct. Later we'll explore other ways to interrogate and visualize
+functions we've represented in our finite element space.
+
+.. doctest::
+
+   >>> fe_approximation[:] = 1  # a function that is 1 everywhere; [:] means "all dofs"
+   >>> plt.subplots(figsize=(6,5))
+   >>> skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True)
+   >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   >>> plt.xlabel('x[0]'); plt.ylabel('x[1]');
+
+.. plot::
+
+   import skfem
+   import matplotlib.pyplot as plt
+   import numpy as np
+   import skfem.visuals.matplotlib
+   
+   mesh = skfem.MeshTri().refined(1)
+   basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
+   fe_approximation = basis_p1.zeros()
+   fe_approximation[:] = 1  # a function that is 1 everywhere; [:] means "all dofs"
+   plt.subplots(figsize=(6,5))
+   skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True)
+   skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   plt.xlabel('x[0]'); plt.ylabel('x[1]');
+
+Now, suppose we want to change this function so it is 0 on the left
+edge. To tell ``skfem`` to make the function zero along those vertical
+line segments on the left edge, we'll call on a very powerful and
+flexible feature of our basis object: ``get_dofs()``.
+
+We can use this method to make ``skfem`` return the indexes to use with
+``fe_approximation`` in order to specify the value of our function in two
+ways: along facets and over entire triangles (``skfem`` calls these
+triangles "cells"/"elements". In this context, "cell"/"element" is purely geometrical
+and should not be confused with the "finite element" which includes a
+concept of polynomial degree.)
+
+.. doctest::
+
+   >>> def is_on_left_edge(x):
+   >>>     return x[0] < 0.1
+   >>> dof_subset_left_edge = basis_p1.get_dofs(facets=is_on_left_edge)
+   >>> fe_approximation[dof_subset_left_edge] = 0
+   >>> plt.subplots(figsize=(6,5))
+   >>> skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
+   >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   >>> plt.xlabel('x[0]'); plt.ylabel('x[1]');
+
+.. plot::
+
+   import skfem
+   import matplotlib.pyplot as plt
+   import numpy as np
+   import skfem.visuals.matplotlib
+   
+   mesh = skfem.MeshTri().refined(1)
+   basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
+   fe_approximation = basis_p1.zeros()
+   fe_approximation[:] = 1  # a function that is 1 everywhere; [:] means "all dofs"
+   def is_on_left_edge(x):
+       return x[0] < 0.1
+   dof_subset_left_edge = basis_p1.get_dofs(facets=is_on_left_edge)
+   fe_approximation[dof_subset_left_edge] = 0
+   plt.subplots(figsize=(6,5))
+   skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
+   skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   plt.xlabel('x[0]'); plt.ylabel('x[1]');
+
+We could make a more complicated function, leaving 0 on that left
+edge, and going to 2 on the top edge. Here we use a lambda function to
+make the code more compact. In general though, lambda functions should
+only be used in trivial circumstances. The verbose naming above is
+more descriptive and readable.
+
+.. doctest::
+
+   >>> dof_subset_right_edge = basis_p1.get_dofs(facets=lambda x: x[1] > 0.9)
+   >>> fe_approximation[dof_subset_right_edge] = 2
+   >>> plt.subplots(figsize=(6,5))
+   >>> skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
+   >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   >>> plt.xlabel('x[0]'); plt.ylabel('x[1]');
+
+.. plot::
+
+   import skfem
+   import matplotlib.pyplot as plt
+   import numpy as np
+   import skfem.visuals.matplotlib
+   
+   mesh = skfem.MeshTri().refined(1)
+   basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
+   fe_approximation = basis_p1.zeros()
+   fe_approximation[:] = 1  # a function that is 1 everywhere; [:] means "all dofs"
+   def is_on_left_edge(x):
+       return x[0] < 0.1
+   dof_subset_left_edge = basis_p1.get_dofs(facets=is_on_left_edge)
+   fe_approximation[dof_subset_left_edge] = 0
+   dof_subset_right_edge = basis_p1.get_dofs(facets=lambda x: x[1] > 0.9)
+   fe_approximation[dof_subset_right_edge] = 2
+   plt.subplots(figsize=(6,5))
+   skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
+   skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   plt.xlabel('x[0]'); plt.ylabel('x[1]');
+
+In a directly analogous manner, we can specify values over entire elements instead of just edges.
+
+.. doctest::
+
+    >>> # reset the function to be 1 everywhere
+    >>> fe_approximation[:] = 1
+    >>> dof_subset_bottom_left = basis_p1.get_dofs(elements=lambda x: np.logical_and(x[0]<.3, x[1]<.3))
+    >>> fe_approximation[dof_subset_bottom_left] = 0
+    >>> plt.subplots(figsize=(6,5))
+    >>> skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
+    >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+    >>> plt.xlabel('x[0]'); plt.ylabel('x[1]');
+
+.. plot::
+
+   import skfem
+   import matplotlib.pyplot as plt
+   import numpy as np
+   import skfem.visuals.matplotlib
+   
+   mesh = skfem.MeshTri().refined(1)
+   basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
+   fe_approximation = basis_p1.zeros()
+   fe_approximation[:] = 1  # a function that is 1 everywhere; [:] means "all dofs"
+   def is_on_left_edge(x):
+       return x[0] < 0.1
+   dof_subset_left_edge = basis_p1.get_dofs(facets=is_on_left_edge)
+   fe_approximation[dof_subset_left_edge] = 0
+   dof_subset_right_edge = basis_p1.get_dofs(facets=lambda x: x[1] > 0.9)
+   fe_approximation[dof_subset_right_edge] = 2
+   # reset the function to be 1 everywhere
+   fe_approximation[:] = 1
+   dof_subset_bottom_left = basis_p1.get_dofs(elements=lambda x: np.logical_and(x[0]<.3, x[1]<.3))
+   fe_approximation[dof_subset_bottom_left] = 0
+   plt.subplots(figsize=(6,5))
+   skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
+   skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   plt.xlabel('x[0]'); plt.ylabel('x[1]');
+
+This is exactly correct. The function is 0 everywhere in the bottom
+left triangle, and goes linearly (because we're in P1) to 1 outside of
+this triangle. Note the continuity between triangles, another
+consequence of using P1 to form our basis set.
+
+To summarize our discussion so far, we've seen how to construct a
+finite element basis set from a mesh and a choice of function over one
+cell of that mesh, in our case P1 (linear polynomials). And we've now
+seen how to create simple functions in that space by specifying the
+value of the function everywhere (``[:]``), along facets
+(``get_dofs(facets=...)``) or over elements (``get_dofs(elements=...)``).
+
+Lets take a closer look at what is happening when we supply a function
+to ``get_dofs()`` by tricking it into plotting the query locations it is
+using. Note the use of lambda here to supply most of the arguments to
+our trial function while still leaving x available as an argument for
+``get_dofs()``.
+
+.. doctest::
+
+   >>> def plot_query_points(x, ax, style, label):
+   >>>     ax.plot(x[0], x[1], style, label=label)
+   >>> plt.subplots(figsize=(5,5))
+   >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   >>> basis_p1.get_dofs(facets=lambda x: plot_query_points(x, plt.gca(), 'or', 'facets'))
+   >>> basis_p1.get_dofs(elements=lambda x: plot_query_points(x, plt.gca(), 'ob', 'elements'))
+   >>> plt.legend()
+
+.. plot::
+
+   import skfem
+   import matplotlib.pyplot as plt
+   import numpy as np
+   import skfem.visuals.matplotlib
+   
+   mesh = skfem.MeshTri().refined(1)
+   basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
+   def plot_query_points(x, ax, style, label):
+       ax.plot(x[0], x[1], style, label=label)
+   plt.subplots(figsize=(5,5))
+   skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   basis_p1.get_dofs(facets=lambda x: plot_query_points(x, plt.gca(), 'or', 'facets'))
+   basis_p1.get_dofs(elements=lambda x: plot_query_points(x, plt.gca(), 'ob', 'elements'))
+   plt.legend()
+
+This plot shows the x coordinates supplied to our test function. If we
+return ``True`` for one of these coordinates, then ``get_dofs()`` will return
+the indexes required by ``fe_approximation`` to force that element or
+facet to a specified value.
+
+The extremely important caveat here is that one should never use ``==``
+when dealing with floating point numbers. Therefore, to find those two
+red dots on the vertical pair of facets in the center, we should write
+as follows. (Later we will show more robust and precise ways of
+labelling facets and elements during mesh construction.)
+
+.. doctest::
+
+   >>> dof_subset_vertical_centerline = basis_p1.get_dofs(facets=lambda x: np.isclose(x[0], 0.5))
+   >>> fe_approximation[:] = 2
+   >>> fe_approximation[dof_subset_vertical_centerline] = 0
+   >>> plt.subplots(figsize=(6,5))
+   >>> skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
+   >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   >>> plt.xlabel('x[0]'); plt.ylabel('x[1]');
+
+.. plot::
+
+   import skfem
+   import matplotlib.pyplot as plt
+   import numpy as np
+   import skfem.visuals.matplotlib
+   
+   mesh = skfem.MeshTri().refined(1)
+   basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
+   fe_approximation = basis_p1.zeros()
+   fe_approximation[:] = 1  # a function that is 1 everywhere; [:] means "all dofs"
+   def is_on_left_edge(x):
+       return x[0] < 0.1
+   dof_subset_left_edge = basis_p1.get_dofs(facets=is_on_left_edge)
+   fe_approximation[dof_subset_left_edge] = 0
+   dof_subset_right_edge = basis_p1.get_dofs(facets=lambda x: x[1] > 0.9)
+   fe_approximation[dof_subset_right_edge] = 2
+   # reset the function to be 1 everywhere
+   fe_approximation[:] = 1
+   dof_subset_bottom_left = basis_p1.get_dofs(elements=lambda x: np.logical_and(x[0]<.3, x[1]<.3))
+   fe_approximation[dof_subset_bottom_left] = 0
+   dof_subset_vertical_centerline = basis_p1.get_dofs(facets=lambda x: np.isclose(x[0], 0.5))
+   fe_approximation[:] = 2
+   fe_approximation[dof_subset_vertical_centerline] = 0
+   plt.subplots(figsize=(6,5))
+   skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
+   skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   plt.xlabel('x[0]'); plt.ylabel('x[1]');
+
+Another way to construct a function in the finite element space is by
+projection. To demonstrate this, we'll use ``f(x) = abs(x[1]-0.5)`` which
+would be a horizontal valley running along the line at ``x[1]=0.5``. We'll
+use an ``skfem`` utility which uses a ``CellBasis`` object to project a Python
+function into the finite element space. The corresponding Python function must
+accept a single argument of point vectors and return an array of
+function values at those points.
+
+.. doctest::
+
+   >>> def f(x):
+   >>>     return 4 * abs(x[1] - 0.5)
+   >>> fe_approximation = basis_p1.project(f)
+   >>> plt.subplots(figsize=(6,5))
+   >>> skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
+   >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   >>> plt.xlabel('x[0]'); plt.ylabel('x[1]');
+
+.. plot::
+
+   import skfem
+   import matplotlib.pyplot as plt
+   import numpy as np
+   import skfem.visuals.matplotlib
+   
+   mesh = skfem.MeshTri().refined(1)
+   basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
+   fe_approximation = basis_p1.zeros()
+   fe_approximation[:] = 1  # a function that is 1 everywhere; [:] means "all dofs"
+   def is_on_left_edge(x):
+       return x[0] < 0.1
+   dof_subset_left_edge = basis_p1.get_dofs(facets=is_on_left_edge)
+   fe_approximation[dof_subset_left_edge] = 0
+   dof_subset_right_edge = basis_p1.get_dofs(facets=lambda x: x[1] > 0.9)
+   fe_approximation[dof_subset_right_edge] = 2
+   # reset the function to be 1 everywhere
+   fe_approximation[:] = 1
+   dof_subset_bottom_left = basis_p1.get_dofs(elements=lambda x: np.logical_and(x[0]<.3, x[1]<.3))
+   fe_approximation[dof_subset_bottom_left] = 0
+   dof_subset_vertical_centerline = basis_p1.get_dofs(facets=lambda x: np.isclose(x[0], 0.5))
+   fe_approximation[:] = 2
+   fe_approximation[dof_subset_vertical_centerline] = 0
+   def f(x):
+       return 4 * abs(x[1] - 0.5)
+   fe_approximation = basis_p1.project(f)
+   plt.subplots(figsize=(6,5))
+   skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
+   skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   plt.xlabel('x[0]'); plt.ylabel('x[1]');
+
+Compare the similarities between this example and the previous one to
+see how there may be more than one way to construct the same function
+in our finite element space. From this point forward, we will refer to
+this process generically as "projecting into the finite element space"
+regardless of which of the methods was actually used to generate the
+projection.
+
+The ``basis_p1`` object and the ``fe_approximation`` array that we've been
+working with are abstract representations of our function in the
+finite element space. Internally, ``skfem`` samples this function at a set
+of locations called "quadrature points". ``skfem`` uses weight sums of
+these samples to compute the integrals it uses to solve PDEs.
+
+These samples at quadrature points are another way to represent the
+functions we have projected into finite element space and it is
+important to understand their relationship with the projections we've
+been constructing. To start this discussion, however, it is important
+to distinguish between "local" coordinates and "global"
+coordinates. In this triangulation we've been working in, the local,
+or reference, triangle is on with vertexes and (0, 0), (1, 0), and (0, 1).
+
+.. doctest::
+
+   >>> plt.subplots(figsize=(5,5))
+   >>> plt.plot([0,1,0,0], [0,0,1,0], 'k')
+   >>> plt.xlabel('x[0] (local coords)'); plt.ylabel('x[1] (local coords)');
+
+.. plot::
+
+   import matplotlib.pyplot as plt
+   plt.subplots(figsize=(5,5))
+   plt.plot([0,1,0,0], [0,0,1,0], 'k')
+   plt.xlabel('x[0] (local coords)'); plt.ylabel('x[1] (local coords)');
+
+Each of the triangles in our mesh can be individually be transformed
+into these coordinates, i.e. for the purposes of integration. The
+quadrature points used are available via the basis object we
+constructed previously, so we can plot their locations on the
+reference triangle.
+
+.. doctest::
+
+   >>> plt.subplots(figsize=(5,5))
+   >>> plt.plot([0,1,0,0], [0,0,1,0], 'k')
+   >>> points, weights = basis_p1.quadrature
+   >>> plt.plot(points[0], points[1], 'or')
+   >>> plt.xlabel('x[0] (local coords)'); plt.ylabel('x[1] (local coords)');
+
+.. plot::
+
+   import skfem
+   import matplotlib.pyplot as plt
+   import numpy as np
+   import skfem.visuals.matplotlib
+   
+   mesh = skfem.MeshTri().refined(1)
+   basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
+   plt.subplots(figsize=(5,5))
+   plt.plot([0,1,0,0], [0,0,1,0], 'k')
+   points, weights = basis_p1.quadrature
+   plt.plot(points[0], points[1], 'or')
+   plt.xlabel('x[0] (local coords)'); plt.ylabel('x[1] (local coords)');
+
+We can get a global visualization of the quadrature points by reverse
+mapping the local coordinates to each of the triangles in our mesh.
+
+.. doctest::
+
+   >>> global_points = basis_p1.mapping.F(points)
+   >>> plt.subplots(figsize=(5,5))
+   >>> plt.plot(global_points[0], global_points[1], 'or')
+   >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   >>> plt.xlabel('x[0]'); plt.ylabel('x[1]');
+
+.. plot::
+
+   import skfem
+   import matplotlib.pyplot as plt
+   import numpy as np
+   import skfem.visuals.matplotlib
+   
+   mesh = skfem.MeshTri().refined(1)
+   basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
+
+   points, weights = basis_p1.quadrature
+   global_points = basis_p1.mapping.F(points)
+   plt.subplots(figsize=(5,5))
+   plt.plot(global_points[0], global_points[1], 'or')
+   skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   plt.xlabel('x[0]'); plt.ylabel('x[1]');
+
+The ``global_points`` array is organized as (coordinate, element_index, quadrature_index):
+
+.. doctest::
+
+   >>> global_points.shape  # 2 dimensional, 8 elements, 3 points/element
+   (2, 8, 3)
+
+To demonstrate how interpolation works, let's annotate each of those
+quadrature points with the values of a (projected) function sampled at
+those locations. To do this, we'll use the ``interpolate`` method of our
+basis object on a function projected into finite element space.
+
+.. doctest::
+
+   >>> def f(x):
+   >>>     return x[0] + x[1]
+   >>> fe_approximation = basis_p1.project(f)
+   >>> interpolation = basis_p1.interpolate(fe_approximation)
+   >>> global_points = basis_p1.mapping.F(points).reshape(2, -1)
+   >>> fig, ax = plt.subplots(1, 2, figsize=(12,6))
+   >>> skfem.visuals.matplotlib.draw(mesh, ax=ax[0])
+   >>> for value, p in zip(interpolation.value.reshape(-1), global_points.T):
+   >>>     ax[0].plot(p[0], p[1], 'or')
+   >>>     ax[0].annotate(f'{value:.2f}', [p[0], p[1]])
+   >>> skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=ax[1], shading='gouraud', colorbar=True)
+   >>> skfem.visuals.matplotlib.draw(mesh, ax=ax[1])
+   >>> ax[1].plot(global_points[0], global_points[1], 'or')
+   >>> plt.xlabel('x[0]'); plt.ylabel('x[1]');
+
+
+.. plot::
+
+   import skfem
+   import matplotlib.pyplot as plt
+   import numpy as np
+   import skfem.visuals.matplotlib
+   
+   mesh = skfem.MeshTri().refined(1)
+   basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
+   points, weights = basis_p1.quadrature
+   global_points = basis_p1.mapping.F(points)
+
+   def f(x):
+       return x[0] + x[1]
+   fe_approximation = basis_p1.project(f)
+   interpolation = basis_p1.interpolate(fe_approximation)
+   global_points = basis_p1.mapping.F(points).reshape(2, -1)
+   fig, ax = plt.subplots(1, 2, figsize=(12,6))
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[0])
+   for value, p in zip(interpolation.value.reshape(-1), global_points.T):
+       ax[0].plot(p[0], p[1], 'or')
+       ax[0].annotate(f'{value:.2f}', [p[0], p[1]])
+   skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=ax[1], shading='gouraud', colorbar=True)
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[1])
+   ax[1].plot(global_points[0], global_points[1], 'or')
+   plt.xlabel('x[0]'); plt.ylabel('x[1]');
+
+The number of quadrature points in an element controls the level of
+accuracy of the integrations. For low degree polynomial basis
+functions, one can supply enough quadrature points for exact
+integration, where the only source of error is the finite machine
+precision of the computer. Using more quadrature points than necessary
+does not further improve accuracy and slightly increases computation
+time, but it can provide a common space to perform computations on
+functions that were projected into different finite element
+spaces. For this reason, it is usually preferred to construct the
+highest order basis set first (in the present consideration that is
+P1) and then derive the lower order basis set from it. This will
+ensure the basis sets share a common set of quadrature points, and
+that there are enough quadrature points to perform exact integration
+of the highest order basis set.
+
+.. doctest::
+
+   >>> basis_p0 = basis_p1.with_element(skfem.ElementTriP0())
+
+The P0 space has functions that are constant over a cell/element in the mesh
+and consequently discontinuous on the facets between cells/elements. It also
+has fewer degrees of freedom than a P1 basis constructed on the same
+mesh. Specifically, the P0 basis will have a degree-of-freedom for
+each cell/element in the mesh.
+
+.. doctest::
+
+   >>> print(f'{basis_p1.zeros().shape[0]} dofs in P1 == {mesh.p.shape[1]} nodes in the mesh')
+   >>> print(f'{basis_p0.zeros().shape[0]} dofs in P0 == {mesh.t.shape[1]} elements in the mesh')
+   9 dofs in P1 == 9 nodes in the mesh
+   8 dofs in P0 == 8 elements in the mesh
+
+Functions can be projected into the P0 space in the same ways that
+were used for P1 projection: ``get_dofs()`` and ``project()``. As the first
+example, we will examine ``get_dofs()`` and compare it to one of the
+previous examples we used in P1: the lower left triangle should be 0
+and 1 everywhere else in the mesh.
diff --git a/docs/index.rst b/docs/index.rst
index 4d964f22..5833db3c 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -27,6 +27,7 @@ Table of contents
 
    self
    gettingstarted
+   extended
    howto
    advanced
    listofexamples

From 9d26ce3c35e11a46096e0154dddfa9668f9dbb5c Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sat, 26 Oct 2024 13:29:25 +0300
Subject: [PATCH 2/6] fix test

---
 docs/extended.rst | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/docs/extended.rst b/docs/extended.rst
index bd106489..244532b6 100644
--- a/docs/extended.rst
+++ b/docs/extended.rst
@@ -1,8 +1,8 @@
 .. _extended:
 
-========================
- Extended documentation
-========================
+===================
+ Extended tutorial
+===================
 
 It is extremely helpful when working with ``skfem`` to understand the
 concepts of quadrature points and degrees-of-freedom, and how these
@@ -325,6 +325,7 @@ our trial function while still leaving x available as an argument for
 
    >>> def plot_query_points(x, ax, style, label):
    >>>     ax.plot(x[0], x[1], style, label=label)
+   >>>     return x[0] * 0
    >>> plt.subplots(figsize=(5,5))
    >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
    >>> basis_p1.get_dofs(facets=lambda x: plot_query_points(x, plt.gca(), 'or', 'facets'))
@@ -342,6 +343,7 @@ our trial function while still leaving x available as an argument for
    basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
    def plot_query_points(x, ax, style, label):
        ax.plot(x[0], x[1], style, label=label)
+       return x[0] * 0
    plt.subplots(figsize=(5,5))
    skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
    basis_p1.get_dofs(facets=lambda x: plot_query_points(x, plt.gca(), 'or', 'facets'))

From e01e612f643e1fc1f3e15ede860bce925b88dcda Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sat, 26 Oct 2024 13:40:37 +0300
Subject: [PATCH 3/6] try avoiding incorrect test failure by sphinx

---
 docs/extended.rst | 14 +++++++++++---
 1 file changed, 11 insertions(+), 3 deletions(-)

diff --git a/docs/extended.rst b/docs/extended.rst
index 244532b6..e6473aa6 100644
--- a/docs/extended.rst
+++ b/docs/extended.rst
@@ -4,6 +4,12 @@
  Extended tutorial
 ===================
 
+.. note::
+
+   This is a tutorial for understanding the classes ``Mesh`` and ``Basis``
+   and how they relate to ``numpy`` arrays used in ``skfem`` to represent
+   finite element functions.
+
 It is extremely helpful when working with ``skfem`` to understand the
 concepts of quadrature points and degrees-of-freedom, and how these
 concepts relate to symbolic math expressions, Python functions, and
@@ -42,7 +48,7 @@ This is easy enough to draw in ``matplotlib``. The first row in the points
 array contains the x locations of each point and the second row
 contains the y locations.
 
-.. doctest::
+.. sourcecode::
 
    >>> plt.plot(mesh.p[0], mesh.p[1], 'ok')
 
@@ -606,7 +612,9 @@ precision of the computer. Using more quadrature points than necessary
 does not further improve accuracy and slightly increases computation
 time, but it can provide a common space to perform computations on
 functions that were projected into different finite element
-spaces. For this reason, it is usually preferred to construct the
+spaces.
+
+For this reason, it is usually preferred to construct the
 highest order basis set first (in the present consideration that is
 P1) and then derive the lower order basis set from it. This will
 ensure the basis sets share a common set of quadrature points, and
@@ -619,7 +627,7 @@ of the highest order basis set.
 
 The P0 space has functions that are constant over a cell/element in the mesh
 and consequently discontinuous on the facets between cells/elements. It also
-has fewer degrees of freedom than a P1 basis constructed on the same
+has fewer degrees-of-freedom than a P1 basis constructed on the same
 mesh. Specifically, the P0 basis will have a degree-of-freedom for
 each cell/element in the mesh.
 

From 41289842a268354b3f7a0b5d38dcb988f87bf328 Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sat, 26 Oct 2024 13:45:34 +0300
Subject: [PATCH 4/6] try skipping doctest

---
 docs/extended.rst | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/docs/extended.rst b/docs/extended.rst
index e6473aa6..7caf7711 100644
--- a/docs/extended.rst
+++ b/docs/extended.rst
@@ -48,9 +48,9 @@ This is easy enough to draw in ``matplotlib``. The first row in the points
 array contains the x locations of each point and the second row
 contains the y locations.
 
-.. sourcecode::
+.. doctest::
 
-   >>> plt.plot(mesh.p[0], mesh.p[1], 'ok')
+   >>> plt.plot(mesh.p[0], mesh.p[1], 'ok')   # doctest: +SKIP
 
 .. plot::
 
@@ -82,8 +82,8 @@ lines to the plot.
        plt.plot(mesh.p[0,[t[2],t[0]]], mesh.p[1,[t[2],t[0]]], 'k')  # from vertex 2 back to 0
 
 This is the most basic triangulation of the unit square. ``skfem`` has
-utilities to refined this mesh into more and smaller triangles as well
-as translate it or map it into different coordinates. It also has
+utilities to refine this mesh into more and smaller triangles as well
+as translate it or map it into different coordinates. It also has some
 basic visualization functions to make drawing the mesh on an ``matplotlib`` axis
 easier.
 

From 0c5d235ad42013009b7425311aca192dc8180e29 Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sat, 26 Oct 2024 13:56:52 +0300
Subject: [PATCH 5/6] remove doctests

---
 docs/extended.rst | 46 +++++++++++++++++++++++-----------------------
 1 file changed, 23 insertions(+), 23 deletions(-)

diff --git a/docs/extended.rst b/docs/extended.rst
index 7caf7711..c280cd22 100644
--- a/docs/extended.rst
+++ b/docs/extended.rst
@@ -48,9 +48,9 @@ This is easy enough to draw in ``matplotlib``. The first row in the points
 array contains the x locations of each point and the second row
 contains the y locations.
 
-.. doctest::
+.. sourcecode::
 
-   >>> plt.plot(mesh.p[0], mesh.p[1], 'ok')   # doctest: +SKIP
+   >>> plt.plot(mesh.p[0], mesh.p[1], 'ok')
 
 .. plot::
 
@@ -63,13 +63,13 @@ together to form triangles. Each column specifes the indexes of the
 points that defined the vertexes of the triangle. We can add these
 lines to the plot.
 
-.. doctest::
+.. sourcecode::
 
-   >>> plt.plot(mesh.p[0], mesh.p[1], 'ok')
+   >>> plt.plot(mesh.p[0], mesh.p[1], 'ok')   # doctest: +SKIP
    >>> for t in mesh.t.T: # transpose to iterate over columns
-   >>>     plt.plot(mesh.p[0,[t[0],t[1]]], mesh.p[1,[t[0],t[1]]], 'k')  # from vertex 0 to 1
-   >>>     plt.plot(mesh.p[0,[t[1],t[2]]], mesh.p[1,[t[1],t[2]]], 'k')  # from vertex 1 to 2
-   >>>     plt.plot(mesh.p[0,[t[2],t[0]]], mesh.p[1,[t[2],t[0]]], 'k')  # from vertex 2 back to 0
+   ...     plt.plot(mesh.p[0,[t[0],t[1]]], mesh.p[1,[t[0],t[1]]], 'k')  # from vertex 0 to 1
+   ...     plt.plot(mesh.p[0,[t[1],t[2]]], mesh.p[1,[t[1],t[2]]], 'k')  # from vertex 1 to 2
+   ...     plt.plot(mesh.p[0,[t[2],t[0]]], mesh.p[1,[t[2],t[0]]], 'k')  # from vertex 2 back to 0
 
 .. plot::
 
@@ -87,7 +87,7 @@ as translate it or map it into different coordinates. It also has some
 basic visualization functions to make drawing the mesh on an ``matplotlib`` axis
 easier.
 
-.. doctest::
+.. sourcecode::
 
    >>> import skfem.visuals.matplotlib
    >>> mesh = skfem.MeshTri().refined(1)
@@ -149,7 +149,7 @@ this will always be equal to the number of nodes in the mesh. However,
 this is in general not true, so to get a ``numpy`` array of the correct
 length and initialized to zeros, we will use our basis object.
 
-.. doctest::
+.. sourcecode::
 
    >>> fe_approximation = basis_p1.zeros()
    
@@ -166,7 +166,7 @@ helper function from ``skfem`` to visualize the functions we
 construct. Later we'll explore other ways to interrogate and visualize
 functions we've represented in our finite element space.
 
-.. doctest::
+.. sourcecode::
 
    >>> fe_approximation[:] = 1  # a function that is 1 everywhere; [:] means "all dofs"
    >>> plt.subplots(figsize=(6,5))
@@ -202,7 +202,7 @@ triangles "cells"/"elements". In this context, "cell"/"element" is purely geomet
 and should not be confused with the "finite element" which includes a
 concept of polynomial degree.)
 
-.. doctest::
+.. sourcecode::
 
    >>> def is_on_left_edge(x):
    >>>     return x[0] < 0.1
@@ -239,7 +239,7 @@ make the code more compact. In general though, lambda functions should
 only be used in trivial circumstances. The verbose naming above is
 more descriptive and readable.
 
-.. doctest::
+.. sourcecode::
 
    >>> dof_subset_right_edge = basis_p1.get_dofs(facets=lambda x: x[1] > 0.9)
    >>> fe_approximation[dof_subset_right_edge] = 2
@@ -272,7 +272,7 @@ more descriptive and readable.
 
 In a directly analogous manner, we can specify values over entire elements instead of just edges.
 
-.. doctest::
+.. sourcecode::
 
     >>> # reset the function to be 1 everywhere
     >>> fe_approximation[:] = 1
@@ -327,7 +327,7 @@ using. Note the use of lambda here to supply most of the arguments to
 our trial function while still leaving x available as an argument for
 ``get_dofs()``.
 
-.. doctest::
+.. sourcecode::
 
    >>> def plot_query_points(x, ax, style, label):
    >>>     ax.plot(x[0], x[1], style, label=label)
@@ -367,7 +367,7 @@ red dots on the vertical pair of facets in the center, we should write
 as follows. (Later we will show more robust and precise ways of
 labelling facets and elements during mesh construction.)
 
-.. doctest::
+.. sourcecode::
 
    >>> dof_subset_vertical_centerline = basis_p1.get_dofs(facets=lambda x: np.isclose(x[0], 0.5))
    >>> fe_approximation[:] = 2
@@ -414,7 +414,7 @@ function into the finite element space. The corresponding Python function must
 accept a single argument of point vectors and return an array of
 function values at those points.
 
-.. doctest::
+.. sourcecode::
 
    >>> def f(x):
    >>>     return 4 * abs(x[1] - 0.5)
@@ -477,7 +477,7 @@ to distinguish between "local" coordinates and "global"
 coordinates. In this triangulation we've been working in, the local,
 or reference, triangle is on with vertexes and (0, 0), (1, 0), and (0, 1).
 
-.. doctest::
+.. sourcecode::
 
    >>> plt.subplots(figsize=(5,5))
    >>> plt.plot([0,1,0,0], [0,0,1,0], 'k')
@@ -496,7 +496,7 @@ quadrature points used are available via the basis object we
 constructed previously, so we can plot their locations on the
 reference triangle.
 
-.. doctest::
+.. sourcecode::
 
    >>> plt.subplots(figsize=(5,5))
    >>> plt.plot([0,1,0,0], [0,0,1,0], 'k')
@@ -522,7 +522,7 @@ reference triangle.
 We can get a global visualization of the quadrature points by reverse
 mapping the local coordinates to each of the triangles in our mesh.
 
-.. doctest::
+.. sourcecode::
 
    >>> global_points = basis_p1.mapping.F(points)
    >>> plt.subplots(figsize=(5,5))
@@ -549,7 +549,7 @@ mapping the local coordinates to each of the triangles in our mesh.
 
 The ``global_points`` array is organized as (coordinate, element_index, quadrature_index):
 
-.. doctest::
+.. sourcecode::
 
    >>> global_points.shape  # 2 dimensional, 8 elements, 3 points/element
    (2, 8, 3)
@@ -559,7 +559,7 @@ quadrature points with the values of a (projected) function sampled at
 those locations. To do this, we'll use the ``interpolate`` method of our
 basis object on a function projected into finite element space.
 
-.. doctest::
+.. sourcecode::
 
    >>> def f(x):
    >>>     return x[0] + x[1]
@@ -621,7 +621,7 @@ ensure the basis sets share a common set of quadrature points, and
 that there are enough quadrature points to perform exact integration
 of the highest order basis set.
 
-.. doctest::
+.. sourcecode::
 
    >>> basis_p0 = basis_p1.with_element(skfem.ElementTriP0())
 
@@ -631,7 +631,7 @@ has fewer degrees-of-freedom than a P1 basis constructed on the same
 mesh. Specifically, the P0 basis will have a degree-of-freedom for
 each cell/element in the mesh.
 
-.. doctest::
+.. sourcecode::
 
    >>> print(f'{basis_p1.zeros().shape[0]} dofs in P1 == {mesh.p.shape[1]} nodes in the mesh')
    >>> print(f'{basis_p0.zeros().shape[0]} dofs in P0 == {mesh.t.shape[1]} elements in the mesh')

From 2d2c4d0cbc208068fba41bfd4b39d409e43b7f7c Mon Sep 17 00:00:00 2001
From: Tom Gustafsson <tom.gustafsson@iki.fi>
Date: Sat, 26 Oct 2024 14:11:06 +0300
Subject: [PATCH 6/6] finish porting

---
 docs/extended.rst | 444 +++++++++++++++++++++++++++++++++++++---------
 1 file changed, 356 insertions(+), 88 deletions(-)

diff --git a/docs/extended.rst b/docs/extended.rst
index c280cd22..f27c8048 100644
--- a/docs/extended.rst
+++ b/docs/extended.rst
@@ -50,7 +50,7 @@ contains the y locations.
 
 .. sourcecode::
 
-   >>> plt.plot(mesh.p[0], mesh.p[1], 'ok')
+   plt.plot(mesh.p[0], mesh.p[1], 'ok')
 
 .. plot::
 
@@ -65,11 +65,11 @@ lines to the plot.
 
 .. sourcecode::
 
-   >>> plt.plot(mesh.p[0], mesh.p[1], 'ok')   # doctest: +SKIP
-   >>> for t in mesh.t.T: # transpose to iterate over columns
-   ...     plt.plot(mesh.p[0,[t[0],t[1]]], mesh.p[1,[t[0],t[1]]], 'k')  # from vertex 0 to 1
-   ...     plt.plot(mesh.p[0,[t[1],t[2]]], mesh.p[1,[t[1],t[2]]], 'k')  # from vertex 1 to 2
-   ...     plt.plot(mesh.p[0,[t[2],t[0]]], mesh.p[1,[t[2],t[0]]], 'k')  # from vertex 2 back to 0
+   plt.plot(mesh.p[0], mesh.p[1], 'ok')
+   for t in mesh.t.T: # transpose to iterate over columns
+       plt.plot(mesh.p[0,[t[0],t[1]]], mesh.p[1,[t[0],t[1]]], 'k')  # from vertex 0 to 1
+       plt.plot(mesh.p[0,[t[1],t[2]]], mesh.p[1,[t[1],t[2]]], 'k')  # from vertex 1 to 2
+       plt.plot(mesh.p[0,[t[2],t[0]]], mesh.p[1,[t[2],t[0]]], 'k')  # from vertex 2 back to 0
 
 .. plot::
 
@@ -89,10 +89,10 @@ easier.
 
 .. sourcecode::
 
-   >>> import skfem.visuals.matplotlib
-   >>> mesh = skfem.MeshTri().refined(1)
-   >>> plt.subplots(figsize=(5,5))
-   >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())  # gca: "get current axis"
+   import skfem.visuals.matplotlib
+   mesh = skfem.MeshTri().refined(1)
+   plt.subplots(figsize=(5,5))
+   skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())  # gca: "get current axis"
 
 .. plot::
 
@@ -151,7 +151,7 @@ length and initialized to zeros, we will use our basis object.
 
 .. sourcecode::
 
-   >>> fe_approximation = basis_p1.zeros()
+   fe_approximation = basis_p1.zeros()
    
 Although this is a simple ``numpy`` array, there are not many things we
 can do with it directly, since out at this level of the code we don't
@@ -168,11 +168,11 @@ functions we've represented in our finite element space.
 
 .. sourcecode::
 
-   >>> fe_approximation[:] = 1  # a function that is 1 everywhere; [:] means "all dofs"
-   >>> plt.subplots(figsize=(6,5))
-   >>> skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True)
-   >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
-   >>> plt.xlabel('x[0]'); plt.ylabel('x[1]');
+   fe_approximation[:] = 1  # a function that is 1 everywhere; [:] means "all dofs"
+   plt.subplots(figsize=(6,5))
+   skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True)
+   skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   plt.xlabel('x[0]'); plt.ylabel('x[1]');
 
 .. plot::
 
@@ -204,14 +204,14 @@ concept of polynomial degree.)
 
 .. sourcecode::
 
-   >>> def is_on_left_edge(x):
-   >>>     return x[0] < 0.1
-   >>> dof_subset_left_edge = basis_p1.get_dofs(facets=is_on_left_edge)
-   >>> fe_approximation[dof_subset_left_edge] = 0
-   >>> plt.subplots(figsize=(6,5))
-   >>> skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
-   >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
-   >>> plt.xlabel('x[0]'); plt.ylabel('x[1]');
+   def is_on_left_edge(x):
+       return x[0] < 0.1
+   dof_subset_left_edge = basis_p1.get_dofs(facets=is_on_left_edge)
+   fe_approximation[dof_subset_left_edge] = 0
+   plt.subplots(figsize=(6,5))
+   skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
+   skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   plt.xlabel('x[0]'); plt.ylabel('x[1]');
 
 .. plot::
 
@@ -241,12 +241,12 @@ more descriptive and readable.
 
 .. sourcecode::
 
-   >>> dof_subset_right_edge = basis_p1.get_dofs(facets=lambda x: x[1] > 0.9)
-   >>> fe_approximation[dof_subset_right_edge] = 2
-   >>> plt.subplots(figsize=(6,5))
-   >>> skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
-   >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
-   >>> plt.xlabel('x[0]'); plt.ylabel('x[1]');
+   dof_subset_right_edge = basis_p1.get_dofs(facets=lambda x: x[1] > 0.9)
+   fe_approximation[dof_subset_right_edge] = 2
+   plt.subplots(figsize=(6,5))
+   skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
+   skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   plt.xlabel('x[0]'); plt.ylabel('x[1]');
 
 .. plot::
 
@@ -274,14 +274,14 @@ In a directly analogous manner, we can specify values over entire elements inste
 
 .. sourcecode::
 
-    >>> # reset the function to be 1 everywhere
-    >>> fe_approximation[:] = 1
-    >>> dof_subset_bottom_left = basis_p1.get_dofs(elements=lambda x: np.logical_and(x[0]<.3, x[1]<.3))
-    >>> fe_approximation[dof_subset_bottom_left] = 0
-    >>> plt.subplots(figsize=(6,5))
-    >>> skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
-    >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
-    >>> plt.xlabel('x[0]'); plt.ylabel('x[1]');
+   # reset the function to be 1 everywhere
+   fe_approximation[:] = 1
+   dof_subset_bottom_left = basis_p1.get_dofs(elements=lambda x: np.logical_and(x[0]<.3, x[1]<.3))
+   fe_approximation[dof_subset_bottom_left] = 0
+   plt.subplots(figsize=(6,5))
+   skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
+   skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   plt.xlabel('x[0]'); plt.ylabel('x[1]');
 
 .. plot::
 
@@ -329,14 +329,14 @@ our trial function while still leaving x available as an argument for
 
 .. sourcecode::
 
-   >>> def plot_query_points(x, ax, style, label):
-   >>>     ax.plot(x[0], x[1], style, label=label)
-   >>>     return x[0] * 0
-   >>> plt.subplots(figsize=(5,5))
-   >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
-   >>> basis_p1.get_dofs(facets=lambda x: plot_query_points(x, plt.gca(), 'or', 'facets'))
-   >>> basis_p1.get_dofs(elements=lambda x: plot_query_points(x, plt.gca(), 'ob', 'elements'))
-   >>> plt.legend()
+   def plot_query_points(x, ax, style, label):
+       ax.plot(x[0], x[1], style, label=label)
+       return x[0] * 0
+   plt.subplots(figsize=(5,5))
+   skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   basis_p1.get_dofs(facets=lambda x: plot_query_points(x, plt.gca(), 'or', 'facets'))
+   basis_p1.get_dofs(elements=lambda x: plot_query_points(x, plt.gca(), 'ob', 'elements'))
+   plt.legend()
 
 .. plot::
 
@@ -369,13 +369,13 @@ labelling facets and elements during mesh construction.)
 
 .. sourcecode::
 
-   >>> dof_subset_vertical_centerline = basis_p1.get_dofs(facets=lambda x: np.isclose(x[0], 0.5))
-   >>> fe_approximation[:] = 2
-   >>> fe_approximation[dof_subset_vertical_centerline] = 0
-   >>> plt.subplots(figsize=(6,5))
-   >>> skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
-   >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
-   >>> plt.xlabel('x[0]'); plt.ylabel('x[1]');
+   dof_subset_vertical_centerline = basis_p1.get_dofs(facets=lambda x: np.isclose(x[0], 0.5))
+   fe_approximation[:] = 2
+   fe_approximation[dof_subset_vertical_centerline] = 0
+   plt.subplots(figsize=(6,5))
+   skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
+   skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   plt.xlabel('x[0]'); plt.ylabel('x[1]');
 
 .. plot::
 
@@ -416,13 +416,13 @@ function values at those points.
 
 .. sourcecode::
 
-   >>> def f(x):
-   >>>     return 4 * abs(x[1] - 0.5)
-   >>> fe_approximation = basis_p1.project(f)
-   >>> plt.subplots(figsize=(6,5))
-   >>> skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
-   >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
-   >>> plt.xlabel('x[0]'); plt.ylabel('x[1]');
+   def f(x):
+       return 4 * abs(x[1] - 0.5)
+   fe_approximation = basis_p1.project(f)
+   plt.subplots(figsize=(6,5))
+   skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=plt.gca(), colorbar=True, shading='gouraud')
+   skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   plt.xlabel('x[0]'); plt.ylabel('x[1]');
 
 .. plot::
 
@@ -479,9 +479,9 @@ or reference, triangle is on with vertexes and (0, 0), (1, 0), and (0, 1).
 
 .. sourcecode::
 
-   >>> plt.subplots(figsize=(5,5))
-   >>> plt.plot([0,1,0,0], [0,0,1,0], 'k')
-   >>> plt.xlabel('x[0] (local coords)'); plt.ylabel('x[1] (local coords)');
+   plt.subplots(figsize=(5,5))
+   plt.plot([0,1,0,0], [0,0,1,0], 'k')
+   plt.xlabel('x[0] (local coords)'); plt.ylabel('x[1] (local coords)');
 
 .. plot::
 
@@ -498,11 +498,11 @@ reference triangle.
 
 .. sourcecode::
 
-   >>> plt.subplots(figsize=(5,5))
-   >>> plt.plot([0,1,0,0], [0,0,1,0], 'k')
-   >>> points, weights = basis_p1.quadrature
-   >>> plt.plot(points[0], points[1], 'or')
-   >>> plt.xlabel('x[0] (local coords)'); plt.ylabel('x[1] (local coords)');
+   plt.subplots(figsize=(5,5))
+   plt.plot([0,1,0,0], [0,0,1,0], 'k')
+   points, weights = basis_p1.quadrature
+   plt.plot(points[0], points[1], 'or')
+   plt.xlabel('x[0] (local coords)'); plt.ylabel('x[1] (local coords)');
 
 .. plot::
 
@@ -524,11 +524,11 @@ mapping the local coordinates to each of the triangles in our mesh.
 
 .. sourcecode::
 
-   >>> global_points = basis_p1.mapping.F(points)
-   >>> plt.subplots(figsize=(5,5))
-   >>> plt.plot(global_points[0], global_points[1], 'or')
-   >>> skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
-   >>> plt.xlabel('x[0]'); plt.ylabel('x[1]');
+   global_points = basis_p1.mapping.F(points)
+   plt.subplots(figsize=(5,5))
+   plt.plot(global_points[0], global_points[1], 'or')
+   skfem.visuals.matplotlib.draw(mesh, ax=plt.gca())
+   plt.xlabel('x[0]'); plt.ylabel('x[1]');
 
 .. plot::
 
@@ -561,20 +561,20 @@ basis object on a function projected into finite element space.
 
 .. sourcecode::
 
-   >>> def f(x):
-   >>>     return x[0] + x[1]
-   >>> fe_approximation = basis_p1.project(f)
-   >>> interpolation = basis_p1.interpolate(fe_approximation)
-   >>> global_points = basis_p1.mapping.F(points).reshape(2, -1)
-   >>> fig, ax = plt.subplots(1, 2, figsize=(12,6))
-   >>> skfem.visuals.matplotlib.draw(mesh, ax=ax[0])
-   >>> for value, p in zip(interpolation.value.reshape(-1), global_points.T):
-   >>>     ax[0].plot(p[0], p[1], 'or')
-   >>>     ax[0].annotate(f'{value:.2f}', [p[0], p[1]])
-   >>> skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=ax[1], shading='gouraud', colorbar=True)
-   >>> skfem.visuals.matplotlib.draw(mesh, ax=ax[1])
-   >>> ax[1].plot(global_points[0], global_points[1], 'or')
-   >>> plt.xlabel('x[0]'); plt.ylabel('x[1]');
+   def f(x):
+       return x[0] + x[1]
+   fe_approximation = basis_p1.project(f)
+   interpolation = basis_p1.interpolate(fe_approximation)
+   global_points = basis_p1.mapping.F(points).reshape(2, -1)
+   fig, ax = plt.subplots(1, 2, figsize=(12,6))
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[0])
+   for value, p in zip(interpolation.value.reshape(-1), global_points.T):
+       ax[0].plot(p[0], p[1], 'or')
+       ax[0].annotate(f'{value:.2f}', [p[0], p[1]])
+   skfem.visuals.matplotlib.plot(basis_p1, fe_approximation, vmin=0, vmax=2, ax=ax[1], shading='gouraud', colorbar=True)
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[1])
+   ax[1].plot(global_points[0], global_points[1], 'or')
+   plt.xlabel('x[0]'); plt.ylabel('x[1]');
 
 
 .. plot::
@@ -623,7 +623,7 @@ of the highest order basis set.
 
 .. sourcecode::
 
-   >>> basis_p0 = basis_p1.with_element(skfem.ElementTriP0())
+   basis_p0 = basis_p1.with_element(skfem.ElementTriP0())
 
 The P0 space has functions that are constant over a cell/element in the mesh
 and consequently discontinuous on the facets between cells/elements. It also
@@ -643,3 +643,271 @@ were used for P1 projection: ``get_dofs()`` and ``project()``. As the first
 example, we will examine ``get_dofs()`` and compare it to one of the
 previous examples we used in P1: the lower left triangle should be 0
 and 1 everywhere else in the mesh.
+
+.. sourcecode::
+
+   # reset the function to be 1 everywhere
+   projection_p1 = basis_p1.zeros()
+   projection_p0 = basis_p0.zeros()
+   projection_p1[:] = 1
+   projection_p0[:] = 1
+   
+   def is_bottom_left(x):
+       return np.logical_and(x[0]<.3, x[1]<.3)
+   
+   p1_bottom_left = basis_p1.get_dofs(elements=is_bottom_left)
+   p0_bottom_left = basis_p0.get_dofs(elements=is_bottom_left)
+   projection_p1[p1_bottom_left] = 0
+   projection_p0[p0_bottom_left] = 0
+   
+   fig, ax = plt.subplots(1, 2, figsize=(12,6))
+   skfem.visuals.matplotlib.plot(basis_p1, projection_p1, vmin=0, vmax=2, ax=ax[0], shading='gouraud', colorbar=True)
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[0])
+   skfem.visuals.matplotlib.plot(basis_p0, projection_p0, vmin=0, vmax=2, ax=ax[1], shading='gouraud', colorbar=True)
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[1])
+   ax[0].set_xlabel('x[0]'); ax[0].set_ylabel('x[1]')
+   ax[1].set_xlabel('x[0]'); ax[1].set_ylabel('x[1]')
+   ax[0].set_title('projection_p1'); ax[1].set_title('projection_p0');
+
+.. plot::
+
+   import skfem
+   import matplotlib.pyplot as plt
+   import numpy as np
+   import skfem.visuals.matplotlib
+   
+   mesh = skfem.MeshTri().refined(1)
+   basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
+   basis_p0 = basis_p1.with_element(skfem.ElementTriP0())
+   # reset the function to be 1 everywhere
+   projection_p1 = basis_p1.zeros()
+   projection_p0 = basis_p0.zeros()
+   projection_p1[:] = 1
+   projection_p0[:] = 1
+   
+   def is_bottom_left(x):
+       return np.logical_and(x[0]<.3, x[1]<.3)
+   
+   p1_bottom_left = basis_p1.get_dofs(elements=is_bottom_left)
+   p0_bottom_left = basis_p0.get_dofs(elements=is_bottom_left)
+   projection_p1[p1_bottom_left] = 0
+   projection_p0[p0_bottom_left] = 0
+   
+   fig, ax = plt.subplots(1, 2, figsize=(12,6))
+   skfem.visuals.matplotlib.plot(basis_p1, projection_p1, vmin=0, vmax=2, ax=ax[0], shading='gouraud', colorbar=True)
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[0])
+   skfem.visuals.matplotlib.plot(basis_p0, projection_p0, vmin=0, vmax=2, ax=ax[1], shading='gouraud', colorbar=True)
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[1])
+   ax[0].set_xlabel('x[0]'); ax[0].set_ylabel('x[1]')
+   ax[1].set_xlabel('x[0]'); ax[1].set_ylabel('x[1]')
+   ax[0].set_title('projection_p1'); ax[1].set_title('projection_p0');
+
+A second example of functions projected into P0 uses ``project()``:
+
+.. sourcecode::
+
+   def f(x):
+       return 2 * abs(x[0] + x[1] - 1)
+   projection_p1 = basis_p1.project(f)
+   projection_p0 = basis_p0.project(f)
+   fig, ax = plt.subplots(1, 2, figsize=(12,6))
+   skfem.visuals.matplotlib.plot(basis_p1, projection_p1, vmin=0, vmax=2, ax=ax[0], shading='gouraud', colorbar=True)
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[0])
+   skfem.visuals.matplotlib.plot(basis_p0, projection_p0, vmin=0, vmax=2, ax=ax[1], shading='gouraud', colorbar=True)
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[1])
+   ax[0].set_xlabel('x[0]'); ax[0].set_ylabel('x[1]')
+   ax[1].set_xlabel('x[0]'); ax[1].set_ylabel('x[1]')
+   ax[0].set_title('projection_p1'); ax[1].set_title('projection_p0');
+
+
+.. plot::
+
+   import skfem
+   import matplotlib.pyplot as plt
+   import numpy as np
+   import skfem.visuals.matplotlib
+   
+   mesh = skfem.MeshTri().refined(1)
+   basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
+   basis_p0 = basis_p1.with_element(skfem.ElementTriP0())
+   def f(x):
+       return 2 * abs(x[0] + x[1] - 1)
+   projection_p1 = basis_p1.project(f)
+   projection_p0 = basis_p0.project(f)
+   fig, ax = plt.subplots(1, 2, figsize=(12,6))
+   skfem.visuals.matplotlib.plot(basis_p1, projection_p1, vmin=0, vmax=2, ax=ax[0], shading='gouraud', colorbar=True)
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[0])
+   skfem.visuals.matplotlib.plot(basis_p0, projection_p0, vmin=0, vmax=2, ax=ax[1], shading='gouraud', colorbar=True)
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[1])
+   ax[0].set_xlabel('x[0]'); ax[0].set_ylabel('x[1]')
+   ax[1].set_xlabel('x[0]'); ax[1].set_ylabel('x[1]')
+   ax[0].set_title('projection_p1'); ax[1].set_title('projection_p0');
+
+It is critical to realize in both of the previous two examples, the
+identical "real" function was projected into each of the P1 and P0
+spaces, and the plots are showing the closest approximation to the
+"real" function available in the respective spaces.
+
+To gain further insight into how well these projections are matching
+our "real" functions, it would be useful to examine these projections
+along a line through the space. For this, we can use the ``probes`` method
+on the basis objects we have constructed. Let's use ``f(x) = 4 * abs(x[0] - 0.5)``
+to illustrate how this works.
+
+.. sourcecode::
+
+   N_query_pts = 100
+   x1_value = 0.25
+   query_pts = np.vstack([
+       np.linspace(0,1,N_query_pts),  # x[0] coordinate values
+       x1_value*np.ones(N_query_pts),  # x[1] coordinate values
+   ])
+   p1_probes = basis_p1.probes(query_pts)
+   p0_probes = basis_p0.probes(query_pts)
+   
+   def f(x):
+       return 4 * abs(x[0] - 0.5)
+   
+   projection_p1 = basis_p1.project(f)
+   projection_p0 = basis_p0.project(f)
+   
+   fig, ax = plt.subplots(2, 2, figsize=(12,12))
+   skfem.visuals.matplotlib.plot(basis_p1, projection_p1, vmin=0, vmax=2, ax=ax[0][0], shading='gouraud')
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[0][0])
+   ax[0][0].plot(query_pts[0], query_pts[1], '--r')
+   skfem.visuals.matplotlib.plot(basis_p0, projection_p0, vmin=0, vmax=2, ax=ax[0][1], shading='gouraud')
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[0][1])
+   ax[0][1].plot(query_pts[0], query_pts[1], '--r')
+   ax[0][0].set_xlabel('x[0]'); ax[0][0].set_ylabel('x[1]')
+   ax[0][1].set_xlabel('x[0]'); ax[0][1].set_ylabel('x[1]')
+   ax[0][0].set_title('projection_p1'); ax[0][1].set_title('projection_p0');
+   
+   ax[1][0].plot(query_pts[0], p1_probes @ projection_p1)
+   ax[1][0].set_xlabel('x[0]'); ax[1][0].set_ylabel('f')
+   ax[1][1].plot(query_pts[0], p0_probes @ projection_p0)
+   ax[1][1].set_xlabel('x[0]'); ax[1][1].set_ylabel('f');
+
+
+.. plot::
+
+   import skfem
+   import matplotlib.pyplot as plt
+   import numpy as np
+   import skfem.visuals.matplotlib
+   
+   mesh = skfem.MeshTri().refined(1)
+   basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
+   basis_p0 = basis_p1.with_element(skfem.ElementTriP0())
+   N_query_pts = 100
+   x1_value = 0.25
+   query_pts = np.vstack([
+       np.linspace(0,1,N_query_pts),  # x[0] coordinate values
+       x1_value*np.ones(N_query_pts),  # x[1] coordinate values
+   ])
+   p1_probes = basis_p1.probes(query_pts)
+   p0_probes = basis_p0.probes(query_pts)
+   
+   def f(x):
+       return 4 * abs(x[0] - 0.5)
+   
+   projection_p1 = basis_p1.project(f)
+   projection_p0 = basis_p0.project(f)
+   
+   fig, ax = plt.subplots(2, 2, figsize=(12,12))
+   skfem.visuals.matplotlib.plot(basis_p1, projection_p1, vmin=0, vmax=2, ax=ax[0][0], shading='gouraud')
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[0][0])
+   ax[0][0].plot(query_pts[0], query_pts[1], '--r')
+   skfem.visuals.matplotlib.plot(basis_p0, projection_p0, vmin=0, vmax=2, ax=ax[0][1], shading='gouraud')
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[0][1])
+   ax[0][1].plot(query_pts[0], query_pts[1], '--r')
+   ax[0][0].set_xlabel('x[0]'); ax[0][0].set_ylabel('x[1]')
+   ax[0][1].set_xlabel('x[0]'); ax[0][1].set_ylabel('x[1]')
+   ax[0][0].set_title('projection_p1'); ax[0][1].set_title('projection_p0');
+   
+   ax[1][0].plot(query_pts[0], p1_probes @ projection_p1)
+   ax[1][0].set_xlabel('x[0]'); ax[1][0].set_ylabel('f')
+   ax[1][1].plot(query_pts[0], p0_probes @ projection_p0)
+   ax[1][1].set_xlabel('x[0]'); ax[1][1].set_ylabel('f');
+
+One more example, using a more refined mesh and a ``f(x) = sin(2 * pi * x[1])``.
+Note that since we are changing the mesh here, we must
+also reconstruct the basis objects.
+
+.. sourcecode::
+
+   mesh = skfem.MeshTri().refined(5)
+   basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
+   basis_p0 = basis_p1.with_element(skfem.ElementTriP0())
+   
+   N_query_pts = 200
+   x0_value = .5
+   query_pts = np.vstack([
+       x0_value * np.ones(N_query_pts),  # x[0] coordinate values
+       np.linspace(0,1,N_query_pts),  # x[1] coordinate values
+   ])
+   p1_probes = basis_p1.probes(query_pts)
+   p0_probes = basis_p0.probes(query_pts)
+   
+   def f(x):
+       return np.sin(2 * np.pi * x[1])
+   
+   projection_p1 = basis_p1.project(f)
+   projection_p0 = basis_p0.project(f)
+   
+   fig, ax = plt.subplots(2, 2, figsize=(12,12))
+   skfem.visuals.matplotlib.plot(basis_p1, projection_p1, vmin=0, vmax=2, ax=ax[0][0], shading='gouraud')
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[0][0])
+   ax[0][0].plot(query_pts[0], query_pts[1], '--r')
+   skfem.visuals.matplotlib.plot(basis_p0, projection_p0, vmin=0, vmax=2, ax=ax[0][1], shading='gouraud')
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[0][1])
+   ax[0][1].plot(query_pts[0], query_pts[1], '--r')
+   ax[0][0].set_xlabel('x[0]'); ax[0][0].set_ylabel('x[1]')
+   ax[0][1].set_xlabel('x[0]'); ax[0][1].set_ylabel('x[1]')
+   ax[0][0].set_title('projection_p1'); ax[0][1].set_title('projection_p0');
+   
+   ax[1][0].plot(query_pts[1], p1_probes @ projection_p1)
+   ax[1][0].set_xlabel('x[1]'); ax[1][0].set_ylabel('f')
+   ax[1][1].plot(query_pts[1], p0_probes @ projection_p0)
+   ax[1][1].set_xlabel('x[1]'); ax[1][1].set_ylabel('f');
+
+.. plot::
+
+   import skfem
+   import matplotlib.pyplot as plt
+   import numpy as np
+   import skfem.visuals.matplotlib
+
+   mesh = skfem.MeshTri().refined(5)
+   basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
+   basis_p0 = basis_p1.with_element(skfem.ElementTriP0())
+   
+   N_query_pts = 200
+   x0_value = .5
+   query_pts = np.vstack([
+       x0_value * np.ones(N_query_pts),  # x[0] coordinate values
+       np.linspace(0,1,N_query_pts),  # x[1] coordinate values
+   ])
+   p1_probes = basis_p1.probes(query_pts)
+   p0_probes = basis_p0.probes(query_pts)
+   
+   def f(x):
+       return np.sin(2 * np.pi * x[1])
+   
+   projection_p1 = basis_p1.project(f)
+   projection_p0 = basis_p0.project(f)
+   
+   fig, ax = plt.subplots(2, 2, figsize=(12,12))
+   skfem.visuals.matplotlib.plot(basis_p1, projection_p1, vmin=0, vmax=2, ax=ax[0][0], shading='gouraud')
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[0][0])
+   ax[0][0].plot(query_pts[0], query_pts[1], '--r')
+   skfem.visuals.matplotlib.plot(basis_p0, projection_p0, vmin=0, vmax=2, ax=ax[0][1], shading='gouraud')
+   skfem.visuals.matplotlib.draw(mesh, ax=ax[0][1])
+   ax[0][1].plot(query_pts[0], query_pts[1], '--r')
+   ax[0][0].set_xlabel('x[0]'); ax[0][0].set_ylabel('x[1]')
+   ax[0][1].set_xlabel('x[0]'); ax[0][1].set_ylabel('x[1]')
+   ax[0][0].set_title('projection_p1'); ax[0][1].set_title('projection_p0');
+   
+   ax[1][0].plot(query_pts[1], p1_probes @ projection_p1)
+   ax[1][0].set_xlabel('x[1]'); ax[1][0].set_ylabel('f')
+   ax[1][1].plot(query_pts[1], p0_probes @ projection_p0)
+   ax[1][1].set_xlabel('x[1]'); ax[1][1].set_ylabel('f');