Merge branch 'feature/addPythonBindings' of ../dune-python into feature/addPythonBindings

dune/python/geometry/CMakeLists.txt

+set(HEADERS
+  multilineargeometry.hh
+  quadraturerules.hh
+  referenceelements.hh
+  type.hh
+)
+
+install(FILES ${HEADERS} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/python/geometry)

dune/python/geometry/multilineargeometry.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#ifndef DUNE_PYTHON_GEOMETRY_MULTILINEARGEOMETRY_HH
+#define DUNE_PYTHON_GEOMETRY_MULTILINEARGEOMETRY_HH
+
+#include <type_traits>
+
+#include <dune/common/visibility.hh>
+
+#include <dune/geometry/type.hh>
+#include <dune/geometry/multilineargeometry.hh>
+
+#include <dune/python/common/typeregistry.hh>
+
+#include <dune/python/pybind11/pybind11.h>
+
+namespace Dune
+{
+
+  namespace Python
+  {
+    // registerMemberFunctions
+    // -----------------------
+
+    template<typename MGeometry>
+    struct RegisterMemberFunctions
+    {
+      RegisterMemberFunctions() {}
+      void operator()(pybind11::class_<MGeometry>& cls)
+      {
+        cls.def("toLocal"                  , &MGeometry::local);
+        cls.def("jacobianInverseTransposed", &MGeometry::jacobianInverseTransposed);
+      }
+    };
+
+
+
+    // doNothing
+    // ---------
+
+    template<typename MGeometry>
+    struct DoNothing
+    {
+      DoNothing() {}
+      void operator()(pybind11::class_<MGeometry>&) {}
+    };
+
+
+
+    // registerMultiLinearGeometryType
+    // -------------------------------
+
+    template<class ctype, int dim, int coorddim>
+    inline auto registerMultiLinearGeometryType(pybind11::module scope)
+    {
+      using pybind11::operator""_a;
+
+      typedef MultiLinearGeometry<ctype, dim, coorddim> MGeometry;
+
+      auto entry = insertClass< MGeometry >( scope, "MultiLinearGeometry_"+std::to_string(dim)+"_"+std::to_string(coorddim),
+          GenerateTypeName("MultiLinearGeometry",MetaType<ctype>(),dim,coorddim),
+          IncludeFiles{"dune/geometry/multilineargeometry.hh"}
+          );
+      auto cls = entry.first;
+
+      if (!entry.second)
+      {
+        cls.def( pybind11::init( [] ( Dune::GeometryType type, pybind11::list corners ) {
+            const std::size_t numCorners = corners.size();
+            std::vector< FieldVector< double, coorddim > > copyCorners( numCorners );
+            for( std::size_t i = 0; i < numCorners; ++i )
+              copyCorners[ i ] = corners[ i ].template cast< FieldVector< double, coorddim > >();
+            return new MGeometry( type, copyCorners );
+          } ), "type"_a, "corners"_a );
+
+        // toLocal and jacobianInverseTransposed call
+        // MatrixHelper::template (xT)rightInvA<m, n> where n has to be >= m (static assert)
+        std::conditional_t<(coorddim >= dim), RegisterMemberFunctions<MGeometry>, DoNothing<MGeometry> >{}(cls);
+
+        cls.def_property_readonly("affine" , [](const MGeometry& geom) { return geom.affine(); });
+        cls.def_property_readonly("type"   , &MGeometry::type);
+        cls.def_property_readonly("corners", &MGeometry::corners);
+        cls.def_property_readonly("center" , &MGeometry::center);
+        cls.def_property_readonly("volume" , &MGeometry::volume);
+        cls.def("corner"                   , &MGeometry::corner);
+        cls.def("integrationElement"       , &MGeometry::integrationElement);
+
+        cls.def("jacobianTransposed",
+            [](const MGeometry& geom, const typename MGeometry::LocalCoordinate& local) {
+              return geom.jacobianTransposed(local);
+            });
+
+        cls.def("toGlobal",
+            [](const MGeometry& geom, const typename MGeometry::LocalCoordinate& local) {
+                return geom.global(local);
+              });
+      }
+#if 0
+      else
+      {
+        scope.def( detail::nameMultiLinearGeometry< dim, coorddim >, [] ( Dune::GeometryType gt, pybind11::list corners ) {
+            std::vector<FieldVector<double, 1>> cornerValues(corners.size());
+            for (unsigned i = 0; i < corners.size(); ++i)
+              cornerValues[i] = corners[i].cast<double>();
+            return MGeometry(gt, cornerValues);
+          } );
+      }
+#endif
+      return cls;
+    }
+
+  } // namespace Python
+
+} // namespace Dune
+
+#endif // ifndef DUNE_PYTHON_GEOMETRY_MULTILINEARGEOMETRY_HH

dune/python/geometry/quadraturerules.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_PYTHON_GEOMETRY_QUADRATURERULES_HH
+#define DUNE_PYTHON_GEOMETRY_QUADRATURERULES_HH
+
+#include <array>
+#include <tuple>
+
+#include <dune/common/visibility.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/geometry/type.hh>
+
+#include <dune/python/common/typeregistry.hh>
+#include <dune/python/pybind11/pybind11.h>
+#include <dune/python/pybind11/numpy.h>
+
+
+namespace Dune
+{
+
+  namespace Python
+  {
+    template <class Rule>
+    auto quadratureToNumpy(const Rule &rule)
+    {
+      pybind11::array_t< double > p(
+          { static_cast< ssize_t >( Rule::d ), static_cast< ssize_t >( rule.size() ) },
+          {
+            static_cast< ssize_t >( reinterpret_cast< const char * >( &rule[ 0 ].position()[ 1 ] ) - reinterpret_cast< const char * >( &rule[ 0 ].position()[ 0 ] ) ),
+            static_cast< ssize_t >( reinterpret_cast< const char * >( &rule[ 1 ].position()[ 0 ] ) - reinterpret_cast< const char * >( &rule[ 0 ].position()[ 0 ] ) )
+          },
+          &rule[ 0 ].position()[ 0 ]
+        );
+      pybind11::array_t< double > w(
+          { static_cast< ssize_t >( rule.size() ) },
+          { static_cast< std::size_t >( reinterpret_cast< const char * >( &rule[ 1 ].weight() ) - reinterpret_cast< const char * >( &rule[ 0 ].weight() ) ) },
+          &rule[ 0 ].weight()
+        );
+
+      return std::make_pair( p, w );
+    }
+
+
+    template <class Rule>
+    auto quadratureToNumpy(pybind11::object self)
+    {
+      const Rule &rule = pybind11::cast< const Rule & >( self );
+      pybind11::array_t< double > p(
+          { static_cast< ssize_t >( Rule::d ), static_cast< ssize_t >( rule.size() ) },
+          {
+            static_cast< ssize_t >( reinterpret_cast< const char * >( &rule[ 0 ].position()[ 1 ] ) - reinterpret_cast< const char * >( &rule[ 0 ].position()[ 0 ] ) ),
+            static_cast< ssize_t >( reinterpret_cast< const char * >( &rule[ 1 ].position()[ 0 ] ) - reinterpret_cast< const char * >( &rule[ 0 ].position()[ 0 ] ) )
+          },
+          &rule[ 0 ].position()[ 0 ],
+          self
+        );
+
+      pybind11::array_t< double > w(
+          { static_cast< ssize_t >( rule.size() ) },
+          { static_cast< std::size_t >( reinterpret_cast< const char * >( &rule[ 1 ].weight() ) - reinterpret_cast< const char * >( &rule[ 0 ].weight() ) ) },
+          &rule[ 0 ].weight(),
+          self
+        );
+
+      return std::make_pair( p, w );
+    }
+    namespace detail
+    {
+
+      // registerQuadraturePoint
+      // -----------------------
+
+      template< class QP >
+      inline void registerQuadraturePoint ( pybind11::object scope, pybind11::class_<QP> cls )
+      {
+        using pybind11::operator""_a;
+
+        typedef typename QP::Vector Vector;
+        typedef typename QP::Field Field;
+
+        cls.def( pybind11::init( [] ( const Vector &x, Field w ) { return new QP( x, w ); } ), "position"_a, "weight"_a );
+
+        cls.def_property_readonly( "position", []( const QP &qp ) -> Vector { return qp.position(); } );
+        cls.def_property_readonly( "weight", &QP::weight );
+
+      }
+
+
+
+      // registerQuadratureRule
+      // ----------------------
+
+      template< class Rule, class... options >
+      inline void registerQuadratureRule ( pybind11::object scope,
+          pybind11::class_<Rule,options...> cls )
+      {
+        cls.def_property_readonly( "order", &Rule::order );
+        cls.def_property_readonly( "type",  &Rule::type );
+
+        cls.def( "get", [] ( pybind11::object self ) {
+            return quadratureToNumpy<Rule>(self);
+            } );
+
+        cls.def( "__iter__", [] ( const Rule &rule ) { return pybind11::make_iterator( rule.begin(), rule.end() ); } );
+      }
+
+    } // namespace detail
+
+
+
+    // registerQuadratureRule
+    // ----------------------
+
+    template< class ctype, int dim >
+    inline auto registerQuadratureRule ( pybind11::object scope )
+    {
+      typedef typename Dune::QuadraturePoint< ctype, dim > QP;
+      auto quadPointCls = insertClass< QP >( scope, "QuadraturePoint",
+            GenerateTypeName("Dune::QuadratePoint",MetaType<ctype>(),dim),
+            IncludeFiles{"dune/python/geometry/quadraturerules.hh"});
+      if (quadPointCls.second)
+        detail::registerQuadraturePoint( scope, quadPointCls.first );
+
+      typedef typename Dune::QuadratureRule< ctype, dim > Rule;
+      auto quadRule = insertClass< Rule >(scope, "QuadratureRule" + std::to_string(dim),
+            GenerateTypeName("Dune::QuadrateRule",MetaType<ctype>(),dim),
+            IncludeFiles{"dune/python/geometry/quadraturerules.hh"});
+      if (quadRule.second)
+        detail::registerQuadratureRule( scope, quadRule.first );
+      return quadRule.first;
+    }
+
+    template< class ctype, int ... dim >
+    inline auto registerQuadratureRule ( pybind11::object scope, std::integer_sequence< int, dim ... > )
+    {
+      return std::make_tuple( registerQuadratureRule< ctype >( scope, std::integral_constant< int, dim >() )... );
+    }
+
+  } // namespace Python
+
+} // namespace Dune
+
+#endif // #ifndef DUNE_PYTHON_GEOMETRY_QUADRATURERULES_HH

dune/python/geometry/referenceelements.hh

+#ifndef DUNE_PYTHON_GEOMETRY_REFERENCEELEMENTS_HH
+#define DUNE_PYTHON_GEOMETRY_REFERENCEELEMENTS_HH
+
+#include <array>
+#include <functional>
+#include <string>
+
+#include <dune/common/visibility.hh>
+
+#include <dune/geometry/referenceelements.hh>
+#include <dune/geometry/type.hh>
+
+#include <dune/python/geometry/quadraturerules.hh>
+#include <dune/python/pybind11/pybind11.h>
+
+namespace Dune
+{
+
+  namespace Python
+  {
+
+    namespace detail
+    {
+
+      // referenceElementSize
+      // --------------------
+
+      template< class RefElement >
+      inline static int referenceElementSize ( const RefElement &refElement, int c )
+      {
+        if( (c < 0) || (c > RefElement::dimension) )
+          throw pybind11::value_error( "Invalid codimension: " + std::to_string( c ) + " (must be in [0, " + std::to_string( RefElement::dimension ) + "])." );
+        return refElement.size( c );
+      }
+
+      template< class RefElement >
+      inline static int referenceElementSize ( const RefElement &refElement, int i, int c, int cc )
+      {
+        const int size = detail::referenceElementSize( refElement, c );
+        if( (i < 0) || (i >= size) )
+          throw pybind11::value_error( "Invalid index: " + std::to_string( i ) + " (must be in [0, " + std::to_string( size ) + "))." );
+        if( (cc < c) || (cc > RefElement::dimension) )
+          throw pybind11::value_error( "Invalid codimension: " + std::to_string( cc ) + " (must be in [" + std::to_string( c ) + ", " + std::to_string( RefElement::dimension ) + "])." );
+        return refElement.size( i, c, cc );
+      }
+
+
+
+      // referenceElementSubEntity
+      // -------------------------
+
+      template< class RefElement >
+      inline static int referenceElementSubEntity ( const RefElement &refElement, int i, int c, int ii, int cc )
+      {
+        const int size = detail::referenceElementSize( refElement, i, c, cc );
+        if( (ii < 0) || (ii >= size) )
+          throw pybind11::value_error( "Invalid index: " + std::to_string( i ) + " (must be in [0, " + std::to_string( size ) + "))." );
+        return refElement.subEntity( i, c, ii, cc );
+      }
+
+
+
+      // referenceElementPosition
+      // ------------------------
+
+      template< class RefElement >
+      inline static auto referenceElementPosition ( const RefElement &refElement, int i, int c )
+      {
+        const int size = detail::referenceElementSize( refElement, c );
+        if( (i < 0) || (i >= size) )
+          throw pybind11::value_error( "Invalid index: " + std::to_string( i ) + " (must be in [0, " + std::to_string( size ) + "))." );
+        return refElement.position( i, c );
+      }
+
+
+      // referenceElementType
+      // --------------------
+
+      template< class RefElement >
+      inline static GeometryType referenceElementType ( const RefElement &refElement, int i, int c )
+      {
+        const int size = detail::referenceElementSize( refElement, c );
+        if( (i < 0) || (i >= size) )
+          throw pybind11::value_error( "Invalid index: " + std::to_string( i ) + " (must be in [0, " + std::to_string( size ) + "))." );
+        return refElement.type( i, c );
+      }
+
+    } // namespace detail
+
+
+    template< class RefElement, class... options >
+    void registerReferenceElements ( pybind11::module module, pybind11::class_< RefElement, options... > cls )
+    {
+      using pybind11::operator""_a;
+
+      pybind11::options opts;
+      opts.disable_function_signatures();
+
+      cls.def_property_readonly( "dimension", [] ( pybind11::object ) { return pybind11::int_( RefElement::dimension ); } );
+
+      cls.def( "size", [] ( const RefElement &self, int c ) {
+          return detail::referenceElementSize( self, c );
+        }, "codim"_a );
+      cls.def( "size", [] ( const RefElement &self, int i, int c, int cc ) {
+          return detail::referenceElementSize( self, i, c, cc );
+         } );
+      cls.def( "size", [] ( const RefElement &self, std::tuple< int, int > e, int cc ) {
+          return detail::referenceElementSize( self, std::get< 0 >( e ), std::get< 1 >( e ), cc );
+        }, "subentity"_a, "codim"_a );
+
+      cls.def( "type", [] ( const RefElement &self, int i, int c ) {
+          return detail::referenceElementType( self, i, c );
+        }, "index"_a, "codim"_a );
+      cls.def( "type", [] ( const RefElement &self, std::tuple< int, int > e ) {
+          return detail::referenceElementType( self, std::get< 0 >( e ), std::get< 1 >( e ) );
+        }, "subentity"_a );
+      cls.def( "types", [] ( const RefElement &self, int c ) {
+          const int size = detail::referenceElementSize( self, c );
+          pybind11::tuple types( size );
+          for( int i = 0; i < size; ++i )
+            types[ i ] = pybind11::cast( self.type( i, c ) );
+          return types;
+        }, "codim"_a,
+        R"doc(
+          get geometry types of subentities
+
+          Args:
+              codim:    codimension of subentities
+
+          Returns:
+              tuple containing geometry type for each subentity
+        )doc" );
+
+      cls.def( "subEntity", [] ( const RefElement &self, int i, int c, int ii, int cc ) {
+          return detail::referenceElementSubEntity( self, i, c, ii, cc );
+        } );
+      cls.def( "subEntity", [] ( const RefElement &self, std::tuple< int, int > e, std::tuple< int, int > ee ) {
+          return detail::referenceElementSubEntity( self, std::get< 0 >( e ), std::get< 1 >( e ), std::get< 0 >( ee ), std::get< 1 >( ee ) );
+        } );
+      cls.def( "subEntities", [] ( const RefElement &self, int i, int c, int cc ) {
+          const int size = detail::referenceElementSize( self, i, c, cc );
+          pybind11::tuple subEntities( size );
+          for( int ii = 0; ii < size; ++ii )
+            subEntities[ ii ] = pybind11::int_( self.subEntity( i, c, ii, cc ) );
+          return subEntities;
+        } );
+      cls.def( "subEntities", [] ( const RefElement &self, std::tuple< int, int > e, int cc ) {
+          const int size = detail::referenceElementSize( self, std::get< 0 >( e ), std::get< 1 >( e ), cc );
+          pybind11::tuple subEntities( size );
+          for( int ii = 0; ii < size; ++ii )
+            subEntities[ ii ] = pybind11::int_( self.subEntity( std::get< 0 >( e ), std::get< 1 >( e ), ii, cc ) );
+          return subEntities;
+        }, "subentity"_a, "codim"_a );
+
+      cls.def( "position", [] ( const RefElement &self, int i, int c ) {
+          return detail::referenceElementPosition( self, i, c );
+        }, "index"_a, "codim"_a );
+      cls.def( "position", [] ( const RefElement &self, std::tuple< int, int > e ) {
+          return detail::referenceElementPosition( self, std::get< 0 >( e ), std::get< 1 >( e ) );
+        }, "subentity"_a );
+      cls.def( "positions", [] ( const RefElement &self, int c ) {
+          const int size = detail::referenceElementSize( self, c );
+          pybind11::tuple positions( size );
+          for( int i = 0; i < size; ++i )
+            positions[ i ] = pybind11::cast( self.position( i, c ) );
+          return positions;
+        }, "codim"_a,
+        R"doc(
+          get barycenters of subentities
+
+          Args:
+              codim:    codimension of subentities
+
+          Returns:
+              tuple containing barycenter for each subentity
+        )doc" );
+
+#if 0
+      // Bug: This property overwrite the method "type"
+      cls.def_property_readonly( "type", [] ( const RefElement &self ) {
+          return self.type();
+        },
+        R"doc(
+          geometry type of reference element
+        )doc" );
+#endif
+      cls.def_property_readonly( "center", [] ( const RefElement &self ) {
+          return self.position( 0, 0 );
+        },
+        R"doc(
+          barycenter of reference domain
+        )doc" );
+      cls.def_property_readonly( "corners", [] ( const RefElement &self ) {
+          const int size = self.size( RefElement::dimension );
+          pybind11::tuple corners( size );
+          for( int i = 0; i < size; ++i )
+            corners[ i ] = pybind11::cast( self.position( i, RefElement::dimension ) );
+          return corners;
+        },
+        R"doc(
+          corners of reference domain
+        )doc" );
+      cls.def_property_readonly( "volume", [] ( const RefElement &self ) {
+          return self.volume();
+        },
+        R"doc(
+          volume of reference domain
+        )doc" );
+
+      if( RefElement::dimension > 0 )
+      {
+        cls.def( "integrationOuterNormal", [] ( const RefElement &self, int i ) {
+            if( (i < 0) || (i >= self.size( 1 )) )
+              throw pybind11::value_error( "Invalid index: " + std::to_string( i ) + " (must be in [0, " + std::to_string( self.size( 1 ) ) + "))." );
+            return self.integrationOuterNormal( i );
+          }, "index"_a );
+        cls.def_property_readonly( "integrationOuterNormals", [] ( const RefElement &self ) {
+            const int size = self.size( 1 );
+            pybind11::tuple integrationOuterNormals( size );
+            for( int i = 0; i < size; ++i )
+              integrationOuterNormals[ i ] = pybind11::cast( self.integrationOuterNormal( i ) );
+            return integrationOuterNormals;
+          } );
+      }
+
+      registerQuadratureRule<typename RefElement::ctype, RefElement::dimension>( module );
+
+      module.def( "general", [] ( const GeometryType &gt ) -> pybind11::object {
+          if( gt.isNone() )
+            return pybind11::none();
+          else
+            return pybind11::cast( Dune::ReferenceElements< typename RefElement::ctype, RefElement::dimension >::general( gt ), pybind11::return_value_policy::reference );
+        } );
+      module.def( "rule", [] (const GeometryType &gt, int order) {
+            return Dune::QuadratureRules< typename RefElement::ctype, RefElement::dimension >::rule( gt, order );
+      }, pybind11::return_value_policy::reference );
+    }
+
+  } // namespace Python
+
+} // namespace Dune
+
+#endif // #ifndef DUNE_PYTHON_GEOMETRY_REFERENCEELEMENTS_HH

dune/python/geometry/type.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_PYTHON_GEOMETRY_TYPE_HH
+#define DUNE_PYTHON_GEOMETRY_TYPE_HH
+
+#include <cassert>
+
+#include <dune/common/exceptions.hh>
+
+#include <dune/geometry/type.hh>
+#include <dune/geometry/typeindex.hh>
+
+#include <dune/python/pybind11/operators.h>
+#include <dune/python/pybind11/pybind11.h>
+
+namespace Dune
+{
+
+  // to_string for GeometryType
+  // --------------------------
+
+  inline static std::string to_string ( const GeometryType &type )
+  {
+    if( type.isNone() )
+      return "none(" + std::to_string( type.dim() ) + ")";
+    switch( type.dim() )
+    {
+    case 0:
+      return "vertex";
+    case 1:
+      return "line";
+    case 2:
+      return (type.isSimplex() ? "triangle" : "quadrilateral");
+    case 3:
+      if( type.isSimplex() )
+        return "tetrahedron";
+      else if( type.isHexahedron() )
+        return "hexahedron";
+      else if( type.isPyramid() )
+        return "pyramid";
+      else if( type.isPrism() )
+        return "prism";
+    default:
+      if( type.isSimplex() )
+        return "simplex(" + std::to_string( type.dim() ) + ")";
+      else if( type.isCube() )
+        return "cube(" + std::to_string( type.dim() ) + ")";
+      else
+        return "general(" + std::to_string( type.id() ) + ", " + std::to_string( type.dim() ) + ")";
+    }
+  }
+
+
+
+  // geometryTypeFromString
+  // ----------------------
+
+  inline static GeometryType geometryTypeFromString ( const std::string &s )
+  {
+    typedef GeometryType (*Constructor) ( const std::vector< std::string > & );
+    static const char *constructorNames[] = {
+        "cube",
+        "general",
+        "hexahedron",
+        "line",
+        "none",
+        "prism",
+        "pyramid",
+        "quadrilateral",
+        "simplex",
+        "tetrahedron",
+        "triangle",
+        "vertex"
+      };
+    static const Constructor constructors[]
+      = {
+          // cube
+          [] ( const std::vector< std::string > &args ) {
+              if( args.size() != 1 )
+                DUNE_THROW( Exception, "GeometryType 'cube' requires integer argument for dimension." );
+              return GeometryTypes::cube( std::stoul( args[ 0 ] ) );
+            },
+          // general
+          [] ( const std::vector< std::string > &args ) {
+              if( args.size() != 2 )
+                DUNE_THROW( Exception, "GeometryType 'general' requires two integer arguments, topologyId and dimension." );
+              const auto id = std::stoul( args[ 0 ] );
+              const auto dim = std::stoul( args[ 1 ] );
+              if( id >= Dune::Impl::numTopologies( dim ) )
+                DUNE_THROW( Exception, "Topology id " << id << " too large for dimension " << dim << "." );
+              return GeometryType( id, dim );
+            },
+          // hexahedron
+          [] ( const std::vector< std::string > &args ) {
+              if( !args.empty() )
+                DUNE_THROW( Exception, "GeometryType 'hexahedron' does not require arguments." );
+              return GeometryTypes::hexahedron;
+            },
+          // line
+          [] ( const std::vector< std::string > &args ) {
+              if( !args.empty() )
+                DUNE_THROW( Exception, "GeometryType 'line' does not require arguments." );
+              return GeometryTypes::line;
+            },
+          // none
+          [] ( const std::vector< std::string > &args ) {
+              if( args.size() != 1 )
+                DUNE_THROW( Exception, "GeometryType 'none' requires integer argument for dimension." );
+              return GeometryTypes::none( std::stoul( args[ 0 ] ) );
+            },
+          // prism
+          [] ( const std::vector< std::string > &args ) {
+              if( !args.empty() )
+                DUNE_THROW( Exception, "GeometryType 'prism' does not require arguments." );
+              return GeometryTypes::prism;
+            },
+          // pyramid
+          [] ( const std::vector< std::string > &args ) {
+              if( !args.empty() )
+                DUNE_THROW( Exception, "GeometryType 'pyramid' does not require arguments." );
+              return GeometryTypes::pyramid;
+            },
+          // quadrilateral
+          [] ( const std::vector< std::string > &args ) {
+              if( !args.empty() )
+                DUNE_THROW( Exception, "GeometryType 'quadrilateral' does not require arguments." );
+              return GeometryTypes::quadrilateral;
+            },
+          // simplex
+          [] ( const std::vector< std::string > &args ) {
+              if( args.size() != 1 )
+                DUNE_THROW( Exception, "GeometryType 'simplex' requires integer argument for dimension." );
+              return GeometryTypes::simplex( std::stoul( args[ 0 ] ) );
+            },
+          // tetrahedron
+          [] ( const std::vector< std::string > &args ) {
+              if( !args.empty() )
+                DUNE_THROW( Exception, "GeometryType 'tetrahedron' does not require arguments." );
+              return GeometryTypes::tetrahedron;
+            },
+          // triangle
+          [] ( const std::vector< std::string > &args ) {
+              if( !args.empty() )
+                DUNE_THROW( Exception, "GeometryType 'triangle' does not require arguments." );
+              return GeometryTypes::triangle;
+            },
+          // vertex
+          [] ( const std::vector< std::string > &args ) {
+              if( !args.empty() )
+                DUNE_THROW( Exception, "GeometryType 'vertex' does not require arguments." );
+              return GeometryTypes::vertex;
+            }
+        };
+    const std::size_t numConstructors = sizeof( constructorNames ) / sizeof( const char * );
+
+    // find constructor index
+    std::size_t n = s.find_first_of( '(' );
+    const std::string cName = s.substr( 0, n );
+    const std::size_t cIdx = std::lower_bound( constructorNames, constructorNames + numConstructors, cName ) - constructorNames;
+    if( (cIdx == numConstructors) || (constructorNames[ cIdx ] != cName) )
+      DUNE_THROW( Exception, "No DUNE geometry type constructor named '" << cName << "'." );
+
+    // obtain argument vector
+    std::vector< std::string > args;
+    if( n != std::string::npos )
+    {
+      while( s[ n ] != ')' )
+      {
+        // skip leading spaces
+        const std::size_t m = s.find_first_not_of( ' ', n+1 );
+        if( m == std::string::npos )
+          DUNE_THROW( Exception, "Invalid argument list." );
+
+        // find end of argument
+        n = s.find_first_of( ",)", m );
+        if( (n == std::string::npos) || (n == m) )
+          DUNE_THROW( Exception, "Invalid argument list." );
+
+        // remove trailing spaces from argument
+        const std::size_t k = s.find_last_not_of( ' ', n-1 );
+        assert( (k != std::string::npos) && (k >= m) );
+
+        args.push_back( s.substr( m, k-m+1 ) );
+      }
+    }
+
+    // call constructor
+    return constructors[ cIdx ]( args );
+  }
+
+
+
+  namespace Python
+  {
+
+    pybind11::class_< GeometryType > registerGeometryType ( pybind11::module scope )
+    {
+      pybind11::class_< GeometryType > cls( scope, "GeometryType" );
+
+      cls.def( pybind11::init( [] ( const std::string &s ) { return new GeometryType( geometryTypeFromString( s ) ); } ) );
+
+      cls.def_property_readonly( "isVertex", [] ( const GeometryType &self ) { return self.isVertex(); } );
+      cls.def_property_readonly( "isLine", [] ( const GeometryType &self ) { return self.isLine(); } );
+      cls.def_property_readonly( "isTriangle", [] ( const GeometryType &self ) { return self.isTriangle(); } );
+      cls.def_property_readonly( "isQuadrilateral", [] ( const GeometryType &self ) { return self.isQuadrilateral(); } );
+      cls.def_property_readonly( "isTetrahedron",[] ( const GeometryType &self ) { return self.isTetrahedron(); } );
+      cls.def_property_readonly( "isPyramid",[] ( const GeometryType &self ) { return self.isPyramid(); } );
+      cls.def_property_readonly( "isPrism", [] ( const GeometryType &self ) { return self.isPrism(); } );
+      cls.def_property_readonly( "isHexahedron", [] ( const GeometryType &self ) { return self.isHexahedron(); } );
+      cls.def_property_readonly( "isSimplex", [] ( const GeometryType &self ) { return self.isSimplex(); } );
+      cls.def_property_readonly( "isCube", [] ( const GeometryType &self ) { return self.isCube(); } );
+      cls.def_property_readonly( "isNone", [] ( const GeometryType &self ) { return self.isNone(); } );
+
+      cls.def( pybind11::self == pybind11::self );
+      cls.def( pybind11::self != pybind11::self );
+      cls.def( "__hash__", [] ( const GeometryType &self ) { return GlobalGeometryTypeIndex::index( self ); } );
+
+      cls.def( "__str__", [] ( const GeometryType &self ) { return to_string( self ); } );
+
+      cls.def_property_readonly( "dim", [] ( const GeometryType &self ) { return self.dim(); } );
+      cls.def_property_readonly( "dimension", [] ( const GeometryType &self ) { return self.dim(); } );
+
+      return cls;
+    }
+
+  } // namespace Python
+
+} // namespace Dune
+
+#endif // ifndef DUNE_PYTHON_GEOMETRY_TYPE_HH