/**
   * \returns 2 for all \p n.
   */
  virtual unsigned int n_second_order_adjacent_vertices (const unsigned int) const override
  { return 2; }

  /**
   * \returns The element-local number of the \f$ v^{th} \f$ vertex

  // adjacent mid-edge shapes. That replaces D by T^-1 diag(T d), which for this T is the sparse
  // update below. alpha = 1/5 is recommended there, and makes the weights strictly positive (QUAD8
  // 1/5 and 4/5, TRI6 1/15 and 1/10) while preserving the partition of unity. T = I elsewhere.
  if (_elem && (_elem->type() == TRI6 || _elem->type() == QUAD8) &&
      get_family() == LAGRANGE && sz == _elem->n_nodes())
    {
      const Real alpha = Real(1)/5;
      // Mid-edge nodes are the trailing indices, and only vertex entries are written, so D(m,m) here
      // is never a value an earlier iteration modified.
      for (const auto m : make_range(_elem->n_vertices(), sz))
        for (const auto v : make_range(_elem->n_second_order_adjacent_vertices(m)))
          {
            const auto vertex = _elem->second_order_adjacent_vertex(m, v);
            D(vertex, vertex) += alpha*D(m,m);
            D(vertex, m)      -= alpha*D(m,m);
          }
    }

unsigned short int Quad8::second_order_adjacent_vertex (const unsigned int n,
                                                        const unsigned int v) const
{
  libmesh_assert_greater_equal (n, this->n_vertices());
  libmesh_assert_less (n, this->n_nodes());
  libmesh_assert_less (v, 2);
  // use the matrix from \p face_quad.C
  return _second_order_adjacent_vertices[n-this->n_vertices()][v];
}