Index: configure.ac
===================================================================
--- configure.ac	(revision 977)
+++ configure.ac	(working copy)
@@ -43,6 +43,7 @@
                 dune/localfunctions/lagrange/pyramidp1/Makefile
                 dune/localfunctions/lagrange/q1/Makefile
                 dune/localfunctions/lagrange/q22d/Makefile
+                dune/localfunctions/lagrange/q23d/Makefile
                 dune/localfunctions/Makefile
                 dune/localfunctions/meta/Makefile
                 dune/localfunctions/meta/power/Makefile
Index: dune/localfunctions/test/Makefile.am
===================================================================
--- dune/localfunctions/test/Makefile.am	(revision 977)
+++ dune/localfunctions/test/Makefile.am	(working copy)
@@ -69,4 +69,7 @@
 TESTS += test-q22d
 test_q22d_SOURCES = test-q22d.cc
 
+TESTS += test-q23d
+test_q23d_SOURCES = test-q23d.cc
+
 include $(top_srcdir)/am/global-rules
Index: dune/localfunctions/test/test-q23d.cc
===================================================================
--- dune/localfunctions/test/test-q23d.cc	(revision 0)
+++ dune/localfunctions/test/test-q23d.cc	(revision 0)
@@ -0,0 +1,59 @@
+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=8 sw=2 sts=2:
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <cstddef>
+#include <iostream>
+#include <ostream>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/geometrytype.hh>
+
+#include <dune/grid/utility/mockgeometry.hh>
+
+#include <dune/localfunctions/lagrange/q23d.hh>
+
+#include "geometries.hh"
+#include "test-fe.hh"
+
+int main(int argc, char** argv) {
+  try {
+    // tolerance for floating-point comparisons
+    static const double eps = 1e-9;
+    // stepsize for numerical differentiation
+    static const double delta = 1e-5;
+
+    int result = 77;
+
+    {
+      std::cout << "== Checking global-valued Q23D elements" << std::endl;
+
+      Dune::GeometryType gt;
+      gt.makeHexahedron();
+
+      typedef TestGeometries<double, 3> TestGeos;
+      static const TestGeos testGeos;
+
+      typedef TestGeos::Geometry Geometry;
+      const Geometry &geo = testGeos.get(gt);
+
+      Dune::Q23DFiniteElementFactory<Geometry, double> feFactory;
+      bool success = testFE(geo, feFactory.make(geo), eps, delta);
+
+      if(success && result != 1)
+        result = 0;
+      else
+        result = 1;
+    }
+
+    return result;
+  }
+  catch (const Dune::Exception& e) {
+    std::cerr << e << std::endl;
+    throw;
+  }
+}
Index: dune/localfunctions/lagrange/q23d.hh
===================================================================
--- dune/localfunctions/lagrange/q23d.hh	(revision 0)
+++ dune/localfunctions/lagrange/q23d.hh	(revision 0)
@@ -0,0 +1,95 @@
+#ifndef DUNE_Q23DLOCALFINITEELEMENT_HH
+#define DUNE_Q23DLOCALFINITEELEMENT_HH
+
+#include <dune/common/geometrytype.hh>
+
+#include <dune/localfunctions/common/localfiniteelementtraits.hh>
+#include <dune/localfunctions/common/localtoglobaladaptors.hh>
+#include "q23d/q23dlocalbasis.hh"
+#include "q23d/q23dlocalcoefficients.hh"
+#include "q23d/q23dlocalinterpolation.hh"
+
+namespace Dune 
+{
+
+  /** \todo Please doc me !
+   */
+  template<class D, class R>
+  class Q23DLocalFiniteElement
+  {
+  public:
+    /** \todo Please doc me !
+     */
+	typedef LocalFiniteElementTraits<Q23DLocalBasis<D,R>,Q23DLocalCoefficients,
+									 Q23DLocalInterpolation<Q23DLocalBasis<D,R> > > Traits;
+
+    /** \todo Please doc me !
+     */
+	Q23DLocalFiniteElement ()
+	{
+	  gt.makeHexahedron();
+	}
+
+    /** \todo Please doc me !
+     */
+	const typename Traits::LocalBasisType& localBasis () const
+	{
+	  return basis;
+	}
+	
+    /** \todo Please doc me !
+     */
+	const typename Traits::LocalCoefficientsType& localCoefficients () const
+	{
+	  return coefficients;
+	}
+	
+    /** \todo Please doc me !
+     */
+	const typename Traits::LocalInterpolationType& localInterpolation () const
+	{
+	  return interpolation;
+	}
+	
+    /** \todo Please doc me !
+     */
+	GeometryType type () const
+	{
+	  return gt;
+	}
+
+  private:
+	Q23DLocalBasis<D,R> basis;
+	Q23DLocalCoefficients coefficients;
+	Q23DLocalInterpolation<Q23DLocalBasis<D,R> > interpolation;
+	GeometryType gt;
+  };
+
+  //! Factory for global-valued Q23D elements
+  /**
+   * \tparam Geometry Type of the geometry.  Used to extract the domain field
+   *                  type.
+   * \tparam RF       Range field type.
+   */
+  template<class Geometry, class RF>
+  class Q23DFiniteElementFactory :
+    public ScalarLocalToGlobalFiniteElementAdaptorFactory<
+      Q23DLocalFiniteElement<typename Geometry::ctype, RF>, Geometry
+      >
+  {
+    typedef Q23DLocalFiniteElement<typename Geometry::ctype, RF> LFE;
+    typedef ScalarLocalToGlobalFiniteElementAdaptorFactory<LFE, Geometry> Base;
+
+    static const LFE lfe;
+
+  public:
+    //! default constructor
+    Q23DFiniteElementFactory() : Base(lfe) {}
+  };
+
+  template<class Geometry, class RF>
+  const typename Q23DFiniteElementFactory<Geometry, RF>::LFE
+  Q23DFiniteElementFactory<Geometry, RF>::lfe;
+}
+
+#endif
Index: dune/localfunctions/lagrange/Makefile.am
===================================================================
--- dune/localfunctions/lagrange/Makefile.am	(revision 977)
+++ dune/localfunctions/lagrange/Makefile.am	(working copy)
@@ -1,11 +1,11 @@
-SUBDIRS = p0 p1 p23d pk2d pk3d prismp1 prismp2 pyramidp1 q1 q22d
+SUBDIRS = p0 p1 p23d pk2d pk3d prismp1 prismp2 pyramidp1 q1 q22d q23d
 
 lagrangedir = $(includedir)/dune/localfunctions/lagrange
 lagrange_HEADERS = emptypoints.hh        lagrangebasis.hh         \
                    equidistantpoints.hh  lagrangecoefficients.hh  \
                    interpolation.hh      lobattopoints.hh         \
                    p0.hh p1.hh p23d.hh p2.hh pk2d.hh pk3d.hh pk.hh pq22d.hh \
-                   prismp1.hh prismp2.hh pyramidp1.hh q1.hh q22d.hh \
+                   prismp1.hh prismp2.hh pyramidp1.hh q1.hh q22d.hh q23d.hh \
                    pqkfactory.hh
 
 check_PROGRAMS = test-all
Index: dune/localfunctions/lagrange/q23d/q23dlocalbasis.hh
===================================================================
--- dune/localfunctions/lagrange/q23d/q23dlocalbasis.hh	(revision 0)
+++ dune/localfunctions/lagrange/q23d/q23dlocalbasis.hh	(revision 0)
@@ -0,0 +1,130 @@
+#ifndef DUNE_Q23DLOCALBASIS_HH
+#define DUNE_Q23DLOCALBASIS_HH
+
+#include <dune/common/fmatrix.hh>
+
+#include <dune/localfunctions/common/localbasis.hh>
+
+namespace Dune 
+{
+/**@ingroup LocalBasisImplementation
+   \brief Lagrange shape functions of order 2 on the reference hexahedron.
+
+   Also known as \f$Q^2\f$.
+
+   \tparam D Type to represent the field in the domain.
+   \tparam R Type to represent the field in the range.
+
+   \nosubgrouping
+*/
+template<class D, class R>
+class Q23DLocalBasis
+{
+public:
+  typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,R,1,Dune::FieldVector<R,1>,
+                             Dune::FieldMatrix<R,1,3> > Traits;
+
+  //! \brief number of shape functions
+  unsigned int size () const
+  {
+	return 27;
+  }
+
+  //! \brief Evaluate all shape functions
+  inline void evaluateFunction (const typename Traits::DomainType& in,
+					 std::vector<typename Traits::RangeType>& out) const
+  { 
+	out.resize(27);
+
+	R x=in[0], y=in[1], z=in[2];
+	R X0=2*x*x-3*x+1, X1=-4*x*x+4*x, X2=2*x*x-x;
+	R Y0=2*y*y-3*y+1, Y1=-4*y*y+4*y, Y2=2*y*y-y;
+	R Z0=2*z*z-3*z+1, Z1=-4*z*z+4*z, Z2=2*z*z-z;
+
+	// numbering: first in descending codim order,
+	// second according to the reference element numbering
+	out[0]  = X0*Y0*Z0;
+	out[14] = X1*Y0*Z0;
+	out[1]  = X2*Y0*Z0;
+	out[12] = X0*Y1*Z0;
+	out[24] = X1*Y1*Z0;
+	out[13] = X2*Y1*Z0;
+	out[2]  = X0*Y2*Z0;
+	out[15] = X1*Y2*Z0;
+	out[3]  = X2*Y2*Z0;
+
+	out[8]  = X0*Y0*Z1;
+	out[22] = X1*Y0*Z1;
+	out[9]  = X2*Y0*Z1;
+	out[20] = X0*Y1*Z1;
+	out[26] = X1*Y1*Z1;
+	out[21] = X2*Y1*Z1;
+	out[10] = X0*Y2*Z1;
+	out[23] = X1*Y2*Z1;
+	out[11] = X2*Y2*Z1;
+
+	out[4]  = X0*Y0*Z2;
+	out[18] = X1*Y0*Z2;
+	out[5]  = X2*Y0*Z2;
+	out[16] = X0*Y1*Z2;
+	out[25] = X1*Y1*Z2;
+	out[17] = X2*Y1*Z2;
+	out[6]  = X0*Y2*Z2;
+	out[19] = X1*Y2*Z2;
+	out[7]  = X2*Y2*Z2;
+  }
+
+  //! \brief Evaluate Jacobian of all shape functions
+  inline void 
+  evaluateJacobian (const typename Traits::DomainType& in,         // position
+			std::vector<typename Traits::JacobianType>& out) const      // return value
+  {  
+	out.resize(27);
+
+	R x=in[0], y=in[1], z=in[2];
+	R X0=2*x*x-3*x+1, X1=-4*x*x+4*x, X2=2*x*x-x;
+	R Y0=2*y*y-3*y+1, Y1=-4*y*y+4*y, Y2=2*y*y-y;
+	R Z0=2*z*z-3*z+1, Z1=-4*z*z+4*z, Z2=2*z*z-z;
+	R DX0=4*x-3, DX1=-8*x+4, DX2=4*x-1;
+	R DY0=4*y-3, DY1=-8*y+4, DY2=4*y-1;
+	R DZ0=4*z-3, DZ1=-8*z+4, DZ2=4*z-1;
+
+	out[0][0][0]  = DX0*Y0*Z0; out[0][0][1]  = X0*DY0*Z0; out[0][0][2]  = X0*Y0*DZ0;
+	out[14][0][0] = DX1*Y0*Z0; out[14][0][1] = X1*DY0*Z0; out[14][0][2] = X1*Y0*DZ0;
+	out[1][0][0]  = DX2*Y0*Z0; out[1][0][1]  = X2*DY0*Z0; out[1][0][2]  = X2*Y0*DZ0;
+	out[12][0][0] = DX0*Y1*Z0; out[12][0][1] = X0*DY1*Z0; out[12][0][2] = X0*Y1*DZ0;
+	out[24][0][0] = DX1*Y1*Z0; out[24][0][1] = X1*DY1*Z0; out[24][0][2] = X1*Y1*DZ0;
+	out[13][0][0] = DX2*Y1*Z0; out[13][0][1] = X2*DY1*Z0; out[13][0][2] = X2*Y1*DZ0;
+	out[2][0][0]  = DX0*Y2*Z0; out[2][0][1]  = X0*DY2*Z0; out[2][0][2]  = X0*Y2*DZ0;
+	out[15][0][0] = DX1*Y2*Z0; out[15][0][1] = X1*DY2*Z0; out[15][0][2] = X1*Y2*DZ0;
+	out[3][0][0]  = DX2*Y2*Z0; out[3][0][1]  = X2*DY2*Z0; out[3][0][2]  = X2*Y2*DZ0;
+
+	out[8][0][0]  = DX0*Y0*Z1; out[8][0][1]  = X0*DY0*Z1; out[8][0][2]  = X0*Y0*DZ1;
+	out[22][0][0] = DX1*Y0*Z1; out[22][0][1] = X1*DY0*Z1; out[22][0][2] = X1*Y0*DZ1;
+	out[9][0][0]  = DX2*Y0*Z1; out[9][0][1]  = X2*DY0*Z1; out[9][0][2]  = X2*Y0*DZ1;
+	out[20][0][0] = DX0*Y1*Z1; out[20][0][1] = X0*DY1*Z1; out[20][0][2] = X0*Y1*DZ1;
+	out[26][0][0] = DX1*Y1*Z1; out[26][0][1] = X1*DY1*Z1; out[26][0][2] = X1*Y1*DZ1;
+	out[21][0][0] = DX2*Y1*Z1; out[21][0][1] = X2*DY1*Z1; out[21][0][2] = X2*Y1*DZ1;
+	out[10][0][0] = DX0*Y2*Z1; out[10][0][1] = X0*DY2*Z1; out[10][0][2] = X0*Y2*DZ1;
+	out[23][0][0] = DX1*Y2*Z1; out[23][0][1] = X1*DY2*Z1; out[23][0][2] = X1*Y2*DZ1;
+	out[11][0][0] = DX2*Y2*Z1; out[11][0][1] = X2*DY2*Z1; out[11][0][2] = X2*Y2*DZ1;
+
+	out[4][0][0]  = DX0*Y0*Z2; out[4][0][1]  = X0*DY0*Z2; out[4][0][2]  = X0*Y0*DZ2;
+	out[18][0][0] = DX1*Y0*Z2; out[18][0][1] = X1*DY0*Z2; out[18][0][2] = X1*Y0*DZ2;
+	out[5][0][0]  = DX2*Y0*Z2; out[5][0][1]  = X2*DY0*Z2; out[5][0][2]  = X2*Y0*DZ2;
+	out[16][0][0] = DX0*Y1*Z2; out[16][0][1] = X0*DY1*Z2; out[16][0][2] = X0*Y1*DZ2;
+	out[25][0][0] = DX1*Y1*Z2; out[25][0][1] = X1*DY1*Z2; out[25][0][2] = X1*Y1*DZ2;
+	out[17][0][0] = DX2*Y1*Z2; out[17][0][1] = X2*DY1*Z2; out[17][0][2] = X2*Y1*DZ2;
+	out[6][0][0]  = DX0*Y2*Z2; out[6][0][1]  = X0*DY2*Z2; out[6][0][2]  = X0*Y2*DZ2;
+	out[19][0][0] = DX1*Y2*Z2; out[19][0][1] = X1*DY2*Z2; out[19][0][2] = X1*Y2*DZ2;
+	out[7][0][0]  = DX2*Y2*Z2; out[7][0][1]  = X2*DY2*Z2; out[7][0][2]  = X2*Y2*DZ2;
+  }
+
+  //! \brief Polynomial order of the shape functions
+  unsigned int order () const
+  {
+	return 2;
+  }
+};
+}
+#endif
Index: dune/localfunctions/lagrange/q23d/q23dlocalcoefficients.hh
===================================================================
--- dune/localfunctions/lagrange/q23d/q23dlocalcoefficients.hh	(revision 0)
+++ dune/localfunctions/lagrange/q23d/q23dlocalcoefficients.hh	(revision 0)
@@ -0,0 +1,55 @@
+#ifndef DUNE_Q23DLOCALCOEFFICIENTS_HH
+#define DUNE_Q23DLOCALCOEFFICIENTS_HH
+
+#include <cstddef>
+#include <iostream>
+#include <vector>
+
+#include <dune/localfunctions/common/localkey.hh>
+
+namespace Dune 
+{
+
+  /**@ingroup LocalLayoutImplementation
+	 \brief Layout map for Q2 elements
+
+	 \nosubgrouping
+     \implements Dune::LocalCoefficientsVirtualImp
+  */
+  class Q23DLocalCoefficients
+  {
+  public:
+	//! \brief Standard constructor
+	Q23DLocalCoefficients () : li(27)
+	{
+	  for (std::size_t i = 0; i < 8; i++)
+		li[i] = LocalKey(i, 3, 0);
+
+	  for (std::size_t i = 8; i < 20; i++)
+		li[i] = LocalKey(i-8, 2, 0);
+
+	  for (std::size_t i = 20; i < 26; i++)
+		li[i] = LocalKey(i-20, 1, 0);
+
+	  li[26] = LocalKey(0, 0, 0);
+	}
+
+	//! number of coefficients
+	std::size_t size () const
+	{
+	  return 27;
+	}
+
+	//! get i'th index
+	const LocalKey& localKey (std::size_t i) const
+	{
+	  return li[i];
+	} 
+
+  private:
+	std::vector<LocalKey> li;
+  };
+
+}
+
+#endif
Index: dune/localfunctions/lagrange/q23d/Makefile.am
===================================================================
--- dune/localfunctions/lagrange/q23d/Makefile.am	(revision 0)
+++ dune/localfunctions/lagrange/q23d/Makefile.am	(revision 0)
@@ -0,0 +1,5 @@
+q23ddir = $(includedir)/dune/localfunctions/lagrange/q23d
+
+q23d_HEADERS = q23dlocalbasis.hh q23dlocalcoefficients.hh q23dlocalinterpolation.hh
+
+include $(top_srcdir)/am/global-rules
Index: dune/localfunctions/lagrange/q23d/q23dlocalinterpolation.hh
===================================================================
--- dune/localfunctions/lagrange/q23d/q23dlocalinterpolation.hh	(revision 0)
+++ dune/localfunctions/lagrange/q23d/q23dlocalinterpolation.hh	(revision 0)
@@ -0,0 +1,54 @@
+#ifndef DUNE_Q23DLOCALINTERPOLATION_HH
+#define DUNE_Q23DLOCALINTERPOLATION_HH
+
+#include <vector>
+
+namespace Dune 
+{
+  template<class LB>
+  class Q23DLocalInterpolation
+  {
+  public:
+
+	//! \brief Local interpolation of a function
+	template<typename F, typename C>
+	void interpolate (const F& f, std::vector<C>& out) const
+	{
+	  typename LB::Traits::DomainType x;
+	  typename LB::Traits::RangeType y;
+	  
+	  out.resize(27);
+	  x[0] = 0.0; x[1] = 0.0; x[2] = 0.0; f.evaluate(x,y); out[0] = y;
+	  x[0] = 0.5; x[1] = 0.0; x[2] = 0.0; f.evaluate(x,y); out[14] = y;
+	  x[0] = 1.0; x[1] = 0.0; x[2] = 0.0; f.evaluate(x,y); out[1] = y;
+	  x[0] = 0.0; x[1] = 0.5; x[2] = 0.0; f.evaluate(x,y); out[12] = y;
+	  x[0] = 0.5; x[1] = 0.5; x[2] = 0.0; f.evaluate(x,y); out[24] = y;
+	  x[0] = 1.0; x[1] = 0.5; x[2] = 0.0; f.evaluate(x,y); out[13] = y;
+	  x[0] = 0.0; x[1] = 1.0; x[2] = 0.0; f.evaluate(x,y); out[2] = y;
+	  x[0] = 0.5; x[1] = 1.0; x[2] = 0.0; f.evaluate(x,y); out[15] = y;
+	  x[0] = 1.0; x[1] = 1.0; x[2] = 0.0; f.evaluate(x,y); out[3] = y;
+
+	  x[0] = 0.0; x[1] = 0.0; x[2] = 0.5; f.evaluate(x,y); out[8] = y;
+	  x[0] = 0.5; x[1] = 0.0; x[2] = 0.5; f.evaluate(x,y); out[22] = y;
+	  x[0] = 1.0; x[1] = 0.0; x[2] = 0.5; f.evaluate(x,y); out[9] = y;
+	  x[0] = 0.0; x[1] = 0.5; x[2] = 0.5; f.evaluate(x,y); out[20] = y;
+	  x[0] = 0.5; x[1] = 0.5; x[2] = 0.5; f.evaluate(x,y); out[26] = y;
+	  x[0] = 1.0; x[1] = 0.5; x[2] = 0.5; f.evaluate(x,y); out[21] = y;
+	  x[0] = 0.0; x[1] = 1.0; x[2] = 0.5; f.evaluate(x,y); out[10] = y;
+	  x[0] = 0.5; x[1] = 1.0; x[2] = 0.5; f.evaluate(x,y); out[23] = y;
+	  x[0] = 1.0; x[1] = 1.0; x[2] = 0.5; f.evaluate(x,y); out[11] = y;
+
+	  x[0] = 0.0; x[1] = 0.0; x[2] = 1.0; f.evaluate(x,y); out[4] = y;
+	  x[0] = 0.5; x[1] = 0.0; x[2] = 1.0; f.evaluate(x,y); out[18] = y;
+	  x[0] = 1.0; x[1] = 0.0; x[2] = 1.0; f.evaluate(x,y); out[5] = y;
+	  x[0] = 0.0; x[1] = 0.5; x[2] = 1.0; f.evaluate(x,y); out[16] = y;
+	  x[0] = 0.5; x[1] = 0.5; x[2] = 1.0; f.evaluate(x,y); out[25] = y;
+	  x[0] = 1.0; x[1] = 0.5; x[2] = 1.0; f.evaluate(x,y); out[17] = y;
+	  x[0] = 0.0; x[1] = 1.0; x[2] = 1.0; f.evaluate(x,y); out[6] = y;
+	  x[0] = 0.5; x[1] = 1.0; x[2] = 1.0; f.evaluate(x,y); out[19] = y;
+	  x[0] = 1.0; x[1] = 1.0; x[2] = 1.0; f.evaluate(x,y); out[7] = y;
+	}
+  };
+}
+
+#endif
