diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 5d8f37f9a173b9a61ccffbb611da8ba909d6fac6..e65636322528cc52687b64cd0acb4155d893ea93 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -7,3 +7,22 @@ include:
before_script:
- . /duneci/bin/duneci-init-job
- duneci-install-module https://gitlab.dune-project.org/core/dune-common.git
+
+debian-11-gcc-9-17-python:
+ image: duneci/debian:11
+ script: duneci-standard-test
+ variables:
+ DUNECI_TOOLCHAIN: gcc-9-17
+ # so we need some variables to build the dune-py module during execution of the first python test:
+ # we need to find the dune mdoule
+ DUNE_CONTROL_PATH: /duneci/modules:$CI_PROJECT_DIR
+ # the position for the dune-py module
+ DUNE_PY_DIR: /duneci/modules/dune-py
+ # during dune-py build this variable is used - do know a way to access
+ # the CMAKE_FLAGS used to build the modules...
+ DUNE_CMAKE_FLAGS: "CC=gcc-9 CXX=g++-9 -DCMAKE_CXX_FLAGS='-std=c++17 -O2 -g -Wall -fdiagnostics-color=always' -DDUNE_ENABLE_PYTHONBINDINGS=ON -DDUNE_MAX_TEST_CORES=4 -DBUILD_SHARED_LIBS=TRUE -DDUNE_PYTHON_INSTALL_LOCATION=none -DCMAKE_POSITION_INDEPENDENT_CODE=TRUE -DCMAKE_DISABLE_FIND_PACKAGE_LATEX=TRUE -DCMAKE_DISABLE_FIND_PACKAGE_Alberta=TRUE -DCMAKE_DISABLE_FIND_PACKAGE_Vc=TRUE -DCMAKE_DISABLE_DOCUMENTATION=TRUE"
+ # cmake flags we use for all dune moudle - important is that build shared libs is set (need some better way of doing this)
+ DUNECI_CMAKE_FLAGS: $DUNE_CMAKE_FLAGS
+ # finally set the python path to all modules
+ PYTHONPATH: /duneci/modules/dune-common/build-cmake/python:/duneci/modules/dune-geometry/build-cmake/python:$CI_PROJECT_DIR/build-cmake/python
+ tags: [duneci]
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 213a4e140ea28b77b37af528996346bccde077df..88060895e173d69068a3d51c989364505ab6f527 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,11 @@
# Master (will become release 2.8)
+- Python bindings have been moved from the `dune-python` module which is now
+ obsolete. To activate Python bindings the CMake flag
+ `DUNE_ENABLE_PYTHONBINDINGS` needs to be turned on (default is off).
+ Furthermore, flags for either shared library or position independent code
+ needs to be used.
+
## Deprecations and removals
- Remove code needed to use reference elements by reference.
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1b015648c09a5e6e7376ca6ded8bce24e4113da3..bb1733698a2559679b7a8305e45ad2a2e84a4fb0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -24,5 +24,11 @@ add_subdirectory("dune")
add_subdirectory("doc")
add_subdirectory("lib")
+# if Python bindings are enabled, include necessary sub directories.
+if( DUNE_ENABLE_PYTHONBINDINGS )
+ add_subdirectory("python")
+ dune_python_install_package(PATH python)
+endif()
+
# finalize the dune project, e.g. generating config.h etc.
finalize_dune_project(GENERATE_CONFIG_H_CMAKE)
diff --git a/dune/CMakeLists.txt b/dune/CMakeLists.txt
index 933ff889d1110658837e180f86b2aaa8717e36cf..3b7d414bda4350bdf24bb5332abc14970eaff750 100644
--- a/dune/CMakeLists.txt
+++ b/dune/CMakeLists.txt
@@ -1 +1,5 @@
add_subdirectory(geometry)
+# if Python bindings are enabled, include necessary sub directories.
+if( DUNE_ENABLE_PYTHONBINDINGS )
+ add_subdirectory("python")
+endif()
diff --git a/dune/python/CMakeLists.txt b/dune/python/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..71d1b23cdddd91e2a4fb072187aa8fc4ace7a7de
--- /dev/null
+++ b/dune/python/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_subdirectory(geometry)
+add_subdirectory(test)
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..3be0d3965c981856f9e320a7a410d906b4642bcf
--- /dev/null
+++ b/dune/python/geometry/referenceelements.hh
@@ -0,0 +1,250 @@
+#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/common/fvecmatregistry.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();
+
+ static const std::size_t dimension = RefElement::dimension;
+ typedef typename RefElement::ctype ctype;
+ registerFieldVecMat<FieldVector<ctype,dimension>>::apply();
+ registerFieldVecMat<FieldMatrix<ctype,dimension,dimension>>::apply();
+
+ 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
diff --git a/dune/python/test/CMakeLists.txt b/dune/python/test/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..2decb23b40e554b112b644ce7b1a9cf4425c82ab
--- /dev/null
+++ b/dune/python/test/CMakeLists.txt
@@ -0,0 +1,12 @@
+dune_python_add_test(NAME pyrefelement
+ COMMAND ${PYTHON_EXECUTABLE} refelement.py
+ WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+ LABELS quick)
+dune_python_add_test(NAME pygeometrytype
+ COMMAND ${PYTHON_EXECUTABLE} geometrytype.py
+ WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+ LABELS quick)
+dune_python_add_test(NAME pytestquad
+ COMMAND ${PYTHON_EXECUTABLE} test_quad.py
+ WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+ LABELS quick)
diff --git a/dune/python/test/geometrytype.py b/dune/python/test/geometrytype.py
new file mode 100644
index 0000000000000000000000000000000000000000..dd43abeb3530ef8505a734eef77b5fc7a2288400
--- /dev/null
+++ b/dune/python/test/geometrytype.py
@@ -0,0 +1,20 @@
+from dune.geometry import *
+
+# make sure t can be reconstructed from str(t)
+for t in (vertex, line, triangle, quadrilateral, tetrahedron, pyramid, prism, hexahedron):
+ assert GeometryType(str(t)) == t
+
+for d in range(10):
+ assert GeometryType(str(simplex(d))) == simplex(d)
+ assert GeometryType(str(cube(d))) == cube(d)
+ assert GeometryType(str(none(d))) == none(d)
+
+# make sure simplices with special names can be constructed by general mechanism
+for d, t in enumerate((vertex, line, triangle, tetrahedron)):
+ assert GeometryType("simplex(" + str(d) + ")") == t
+ assert GeometryType("general( 0, " + str(d) + ")") == t
+
+# make sure cube with special names can be constructed by general mechanism
+for d, t in enumerate((vertex, line, quadrilateral, hexahedron)):
+ assert GeometryType("cube(" + str(d) + ")") == t
+ assert GeometryType("general( " + str(2**d-1) + ", " + str(d) + ")") == t
diff --git a/dune/python/test/refelement.py b/dune/python/test/refelement.py
new file mode 100644
index 0000000000000000000000000000000000000000..74ca8d2ac5fac7618ec62b33ff5f3e1b13ce2cd6
--- /dev/null
+++ b/dune/python/test/refelement.py
@@ -0,0 +1,40 @@
+from dune.geometry import referenceElement
+from dune.geometry import vertex, line, triangle, quadrilateral, tetrahedron, pyramid, prism, hexahedron, none
+
+def test(r):
+ for codim in range(r.dimension+1):
+ types = r.types(codim)
+ if len(types) != r.size(codim):
+ raise Exception("types tuple has wrong size")
+ for i in range(len(types)):
+ if types[i] != r.type(i, codim):
+ raise Exception("types tuple has wrong content")
+
+ for codim in range(r.dimension+1):
+ positions = r.positions(codim)
+ if len(positions) != r.size(codim):
+ raise Exception("positions tuple has wrong size")
+ for i in range(len(positions)):
+ if positions[i] != r.position(i, codim):
+ raise Exception("positions tuple has wrong content")
+
+ if r.dimension > 0:
+ normals = r.integrationOuterNormals
+ if len(normals) != r.size(1):
+ raise Exception("integrationOuterNormals has wrong size")
+ for i in range(len(normals)):
+ if normals[i] != r.integrationOuterNormal(i):
+ raise Exception("integrationOuterNormals has wrong content")
+
+test(referenceElement(vertex))
+test(referenceElement(line))
+test(referenceElement(triangle))
+test(referenceElement(quadrilateral))
+test(referenceElement(tetrahedron))
+test(referenceElement(pyramid))
+test(referenceElement(prism))
+test(referenceElement(hexahedron))
+
+for dim in range(4):
+ if referenceElement(none(dim)) is not None:
+ raise Exception("got reference element for geometry type none")
diff --git a/dune/python/test/test_quad.py b/dune/python/test/test_quad.py
new file mode 100644
index 0000000000000000000000000000000000000000..466a8c13d96a2f719fbf52fae3bd6ce795d9432c
--- /dev/null
+++ b/dune/python/test/test_quad.py
@@ -0,0 +1,44 @@
+import time, math, numpy
+import dune.geometry as geo
+
+def monomial(p):
+ def function(point):
+ return sum( x**p for x in point)
+ return function
+
+result = {3: # integral for sum_i x_i^p over reference element
+ {geo.line: 1./4.,
+ geo.triangle: 0.1,
+ geo.quadrilateral: 1./2.,
+ geo.tetrahedron: 1./40.,
+ geo.pyramid: None,
+ geo.prism: None,
+ geo.hexahedron: 3./4.,
+ },
+ 4:
+ {geo.line: 1./5.,
+ geo.triangle: 1./15.,
+ geo.quadrilateral: 2./5.,
+ geo.tetrahedron: 1./70.,
+ geo.pyramid: None,
+ geo.prism: None,
+ geo.hexahedron: 3./5.,
+ },
+ }
+
+
+
+for order in [3,4]:
+ rules = geo.quadratureRules(order)
+ p = monomial(order)
+ for t in (geo.line, geo.triangle, geo.quadrilateral,
+ geo.tetrahedron, geo.pyramid, geo.prism, geo.hexahedron):
+ value1 = 0
+ for q in rules(t):
+ value1 += p(q.position)*q.weight
+ hatxs, hatws = rules(t).get()
+ value2 = numpy.sum(p(hatxs) * hatws, axis=-1)
+ # print(order,t,value2)
+ assert abs(value1-value2)<1e-14
+ if result[order][t] is not None:
+ assert abs(result[order][t] - value1)<1e-12
diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..bf6f08c0410bb4c3fcac06cfadfb8be2eeb46376
--- /dev/null
+++ b/python/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_subdirectory(dune)
+
+configure_file(setup.py.in setup.py)
diff --git a/python/dune/CMakeLists.txt b/python/dune/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a57867b275325a42830280415ddb98983025e42e
--- /dev/null
+++ b/python/dune/CMakeLists.txt
@@ -0,0 +1,4 @@
+add_subdirectory(geometry)
+add_python_targets(dune
+ __init__
+)
diff --git a/python/dune/__init__.py b/python/dune/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..de40ea7ca058e07a399acf61529d418c07eeee61
--- /dev/null
+++ b/python/dune/__init__.py
@@ -0,0 +1 @@
+__import__('pkg_resources').declare_namespace(__name__)
diff --git a/python/dune/geometry/CMakeLists.txt b/python/dune/geometry/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..87cdacc35974670fd9db21b2d36bb09db17bd14e
--- /dev/null
+++ b/python/dune/geometry/CMakeLists.txt
@@ -0,0 +1,7 @@
+add_python_targets(geometry
+ __init__
+ _referenceelements
+ quadpy
+)
+dune_add_pybind11_module(NAME _geometry)
+set_property(TARGET _geometry PROPERTY LINK_LIBRARIES dunecommon dunegeometry)
diff --git a/python/dune/geometry/__init__.py b/python/dune/geometry/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..857d92734bbf56d26b6b6159244301794d3fdf4a
--- /dev/null
+++ b/python/dune/geometry/__init__.py
@@ -0,0 +1,30 @@
+from ._geometry import *
+from ._referenceelements import *
+import numpy
+
+def _duneIntegrate(self,entity,f):
+ points,weights = self.get()
+ try:
+ ie = entity.geometry.integrationElement
+ except AttributeError:
+ ie = geometry.integrationElement
+ return numpy.sum(f(entity,points)*ie(points)*weights,axis=-1)
+
+_duneQuadratureRules = {}
+def quadratureRule(geometryType, order):
+ try:
+ geometryType = geometryType.type
+ except AttributeError:
+ pass
+ try:
+ rule = _duneQuadratureRules[(geometryType,order)]
+ except KeyError:
+ rule = module(geometryType.dim).rule(geometryType,order)
+ setattr(rule.__class__,"apply",_duneIntegrate)
+ _duneQuadratureRules[(geometryType,order)] = rule
+ return rule
+def quadratureRules(order):
+ return lambda entity: quadratureRule(entity,order)
+
+def integrate(rules,entity,f):
+ return rules(entity).apply(entity,f)
diff --git a/python/dune/geometry/_geometry.cc b/python/dune/geometry/_geometry.cc
new file mode 100644
index 0000000000000000000000000000000000000000..defba1f4b21dbc5856a5597b76af675e9bd56806
--- /dev/null
+++ b/python/dune/geometry/_geometry.cc
@@ -0,0 +1,25 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#include <dune/python/geometry/type.hh>
+#include <dune/python/pybind11/pybind11.h>
+
+PYBIND11_MODULE( _geometry, module )
+{
+ Dune::Python::registerGeometryType( module );
+
+ // register geometry type contuctors
+ module.def( "simplex", [] ( int dim ) { return Dune::GeometryTypes::simplex( dim ); } );
+ module.def( "cube", [] ( int dim ) { return Dune::GeometryTypes::cube( dim ); } );
+ module.def( "none", [] ( int dim ) { return Dune::GeometryTypes::none( dim ); } );
+
+ // register predefined geometry types
+ module.attr( "vertex" ) = Dune::GeometryTypes::vertex;
+ module.attr( "line" ) = Dune::GeometryTypes::line;
+ module.attr( "triangle" ) = Dune::GeometryTypes::triangle;
+ module.attr( "quadrilateral" ) = Dune::GeometryTypes::quadrilateral;
+ module.attr( "tetrahedron" ) = Dune::GeometryTypes::tetrahedron;
+ module.attr( "pyramid" ) = Dune::GeometryTypes::pyramid;
+ module.attr( "prism" ) = Dune::GeometryTypes::prism;
+ module.attr( "hexahedron" ) = Dune::GeometryTypes::hexahedron;
+}
diff --git a/python/dune/geometry/_referenceelements.py b/python/dune/geometry/_referenceelements.py
new file mode 100644
index 0000000000000000000000000000000000000000..184706ba8bbed4d709ec7d87b07f169c8f9897bc
--- /dev/null
+++ b/python/dune/geometry/_referenceelements.py
@@ -0,0 +1,22 @@
+from ..generator.generator import SimpleGenerator
+from dune.common.hashit import hashIt
+def module(dim):
+ typeName = "Dune::Geo::ReferenceElement<Dune::Geo::ReferenceElementImplementation<double," + str(dim) + "> >"
+ includes = ["dune/python/geometry/referenceelements.hh"]
+ typeHash = "referenceelements_" + hashIt(typeName)
+ generator = SimpleGenerator("ReferenceElements", "Dune::Python")
+ m = generator.load(includes, typeName, typeHash)
+ return m
+
+_duneReferenceElements = {}
+def referenceElement(geometryType):
+ try:
+ geometryType = geometryType.type
+ except:
+ pass
+ try:
+ ref = _duneReferenceElements[geometryType]
+ except KeyError:
+ ref = module(geometryType.dim).general(geometryType)
+ _duneReferenceElements[geometryType] = ref
+ return ref
diff --git a/python/dune/geometry/quadpy.py b/python/dune/geometry/quadpy.py
new file mode 100644
index 0000000000000000000000000000000000000000..8ba8c3c00f81ef2cc07afb0241457589f591a54b
--- /dev/null
+++ b/python/dune/geometry/quadpy.py
@@ -0,0 +1,111 @@
+import logging
+
+logger = logging.getLogger(__name__)
+
+try:
+ import quadpy as qp
+ import numpy
+ class QPQuadPoint:
+ def __init__(self,p,w):
+ self.p_ = p
+ self.w_ = w
+ @property
+ def position(self):
+ return self.p_
+ @property
+ def weight(self):
+ return self.w_
+ _cache = {}
+ class QPQuadPyQuadrature:
+ def __init__(self,quad,order,method,vertices,transform):
+ self.quad_ = quad
+ self.order_ = order
+ self.method_ = method
+ self.vertices_ = numpy.array(vertices)
+ # here an error will occur if the method is invalid - how to catch?
+ self.quadrature_ = getattr(quad,method)
+ try:
+ self.points_ = transform.transform(self.quadrature_.points.T, self.vertices_)
+ except ValueError:
+ self.points_ = transform.transform(self.quadrature_.points.T, self.vertices_.T).T
+ try:
+ self.weights_ = transform.get_detJ(self.quadrature_.points.T, self.vertices_)*self.quadrature_.weights
+ except AttributeError:
+ self.weights_ = transform.get_vol(self.vertices_)*self.quadrature_.weights
+ self.quadPoints_ = [ QPQuadPoint(p,w) for p,w in zip(self.points_,self.weights_) ]
+ self.points_ = self.points_.transpose().copy()
+ def get(self):
+ return self.points_, self.weights_
+ def apply(self,entity,f):
+ ie = entity.geometry.integrationElement
+ f_e = lambda x: f(entity,x)*ie(x)
+ return self.quad_.integrate(f_e,self.vertices_,self.quadrature_)
+ def __iter__(self):
+ return self.quadPoints_.__iter__()
+ @property
+ def order(self):
+ return self.quadrature_.degree
+ def name(self):
+ return self.method_
+ def rule(gt,quadDescription):
+ try:
+ gt = gt.type
+ except AttributeError:
+ pass
+ try:
+ quad = _cache[(gt,quadDescription)]
+ return quad
+ except KeyError:
+ pass
+ order, method = quadDescription
+ dim = gt.dim
+ if gt.isLine:
+ vertices = [0,1]
+ quad = qp.line_segment
+ from quadpy.ncube import transform
+ elif gt.isTriangle:
+ vertices = [[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]
+ quad = qp.triangle
+ from quadpy.nsimplex import transform
+ elif gt.isQuadrilateral:
+ vertices = qp.quadrilateral.rectangle_points([0.0, 1.0], [0.0, 1.0])
+ quad = qp.quadrilateral
+ from quadpy.ncube import transform
+ elif gt.isTetrahedron:
+ vertices = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
+ quad = qp.tetrahedron
+ from quadpy.nsimplex import transform
+ elif gt.isHexahedron:
+ vertices = qp.hexahedron.cube_points([0.0, 1.0], [0.0, 1.0], [0.0, 1.0])
+ quad = qp.hexahedron
+ from quadpy.ncube import transform
+ elif gt.isPyramid:
+ vertices = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 1.0, 0.0] [0.0, 0.0, 1.0]]
+ quad = qp.pyramid
+ raise ValueError("prism quadratures not yet fully supported")
+ elif gt.isPrism:
+ vertices = [
+ [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0],
+ [0.0, 0.0, 1.0], [1.0, 0.0, 1.0], [0.0, 1.0, 1.0]
+ ]
+ quad = qp.wedge
+ raise ValueError("prism quadratures not yet fully supported")
+ else:
+ raise ValueError("no quadpy quadrature available for the geometryType " + str(gt))
+ if not method and not order:
+ return quad
+ ret = QPQuadPyQuadrature(quad,order,method,vertices,transform)
+ _cache[(gt,quadDescription)] = ret
+ return ret
+
+ def rules(methods):
+ def r(entity):
+ try:
+ return rule(entity.type, methods[entity.type])
+ except AttributeError:
+ return rule(entity, methods[entity])
+ return r
+
+except ImportError as e:
+ logger.warning('Unable to import quadpy: ' + " ".join(str(e).splitlines()))
+ raise ImportError("Unable to import quadpy module")
diff --git a/python/setup.py.in b/python/setup.py.in
new file mode 100644
index 0000000000000000000000000000000000000000..769795bf233e7b9811b71d349535dd66264afb03
--- /dev/null
+++ b/python/setup.py.in
@@ -0,0 +1,12 @@
+from setuptools import setup, find_packages
+
+setup(name="dune.geometry",
+ namespace_packages=['dune'],
+ description="Python lib for dune",
+ version="${DUNE_GEOMETRY_VERSION}",
+ author="Andreas Dedner and Martin Nolte",
+ packages = find_packages(),
+ zip_safe = 0,
+ package_data = {'': ['*.so']},
+ install_requires = ['portalocker']
+ )