diff --git a/dune/python/geometry/CMakeLists.txt b/dune/python/geometry/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f44ebfd0ab464a193168d7172d6699708750ac05
--- /dev/null
+++ b/dune/python/geometry/CMakeLists.txt
@@ -0,0 +1,8 @@
+set(HEADERS
+ multilineargeometry.hh
+ quadraturerules.hh
+ referenceelements.hh
+ type.hh
+)
+
+install(FILES ${HEADERS} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/python/geometry)
diff --git a/dune/python/geometry/multilineargeometry.hh b/dune/python/geometry/multilineargeometry.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d42f3cef5a0420a3c5d1e2035c85bb49ea59862e
--- /dev/null
+++ b/dune/python/geometry/multilineargeometry.hh
@@ -0,0 +1,117 @@
+// -*- 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
diff --git a/dune/python/geometry/quadraturerules.hh b/dune/python/geometry/quadraturerules.hh
new file mode 100644
index 0000000000000000000000000000000000000000..da4c4c7b3e5a7f6b3f940d0051a5c2977268f633
--- /dev/null
+++ b/dune/python/geometry/quadraturerules.hh
@@ -0,0 +1,144 @@
+// -*- 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
diff --git a/dune/python/geometry/referenceelements.hh b/dune/python/geometry/referenceelements.hh
new file mode 100644
index 0000000000000000000000000000000000000000..91339dfe250b478870b491faeb87847840869876
--- /dev/null
+++ b/dune/python/geometry/referenceelements.hh
@@ -0,0 +1,243 @@
+#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 > ) -> 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 >, 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
diff --git a/dune/python/geometry/type.hh b/dune/python/geometry/type.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8398136b6f97b161952393e47c8c4fdb7cb78639
--- /dev/null
+++ b/dune/python/geometry/type.hh
@@ -0,0 +1,230 @@
+// -*- 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