1 #ifndef DUNE_PYTHON_GRID_NUMPY_HH 2 #define DUNE_PYTHON_GRID_NUMPY_HH 3 4 #include <cassert> 5 #include <cstddef> 6 7 #include <algorithm> 8 #include <array> 9 #include <map> 10 #include <vector> 11 12 #include <dune/common/ftraits.hh> 13 14 #include <dune/geometry/dimension.hh> 15 #include <dune/geometry/type.hh> 16 #include <dune/geometry/virtualrefinement.hh> 17 #include <dune/geometry/referenceelements.hh> 18 19 #include <dune/grid/common/mcmgmapper.hh> 20 #include <dune/grid/common/rangegenerators.hh> 21 22 #include <dune/python/common/getdimension.hh> 23 #include <dune/python/grid/object.hh> 24 #include <dune/python/pybind11/numpy.h> 25 #include <dune/python/pybind11/pybind11.h> 26 27 namespace Dune 28 { 29 30 namespace Python 31 { 32 33 // External Forward Declarations 34 // ----------------------------- 35 36 template< class GridFunction > 37 struct GridFunctionTraits; 38 39 40 41 // coordinates 42 // ----------- 43 44 template< class GridView, class Mapper > coordinates(const GridView & gridView,const Mapper & mapper)45 inline static pybind11::array_t< typename GridView::ctype > coordinates ( const GridView &gridView, const Mapper &mapper ) 46 { 47 typedef typename GridView::ctype ctype; 48 49 const std::vector< std::size_t > shape{ static_cast< std::size_t >( mapper.size() ), static_cast< std::size_t >( GridView::dimensionworld ) }; 50 const std::vector< std::size_t > stride{ GridView::dimensionworld * sizeof( ctype ), sizeof( ctype ) }; 51 pybind11::array_t< ctype > coords( pybind11::buffer_info( nullptr, sizeof( ctype ), pybind11::format_descriptor< ctype >::value, 2, shape, stride ) ); 52 53 pybind11::buffer_info info = coords.request( true ); 54 for( const auto &vertex : vertices( gridView, Partitions::all ) ) 55 { 56 typename Mapper::Index index; 57 if( !mapper.contains( vertex, index ) ) 58 continue; 59 60 const auto x = vertex.geometry().center(); 61 std::copy( x.begin(), x.end(), static_cast< ctype * >( info.ptr ) + GridView::dimensionworld * index ); 62 } 63 64 return coords; 65 } 66 67 template< class GridView > coordinates(const GridView & gridView)68 inline static pybind11::array_t< typename GridView::ctype > coordinates ( const GridView &gridView ) 69 { 70 MultipleCodimMultipleGeomTypeMapper< GridView > mapper( gridView, mcmgVertexLayout() ); 71 return coordinates( gridView, mapper ); 72 } 73 74 75 76 // flatCopy 77 // -------- 78 79 template< class In, class T, class = std::enable_if_t< std::is_convertible< In, T >::value > > flatCopy(const In & in,T * out)80 T *flatCopy ( const In &in, T *out ) 81 { 82 *out = in; 83 return ++out; 84 } 85 86 template< class In, class T, class = decltype( std::declval< const In & >().begin() ), class = std::enable_if_t< !std::is_convertible< In, T >::value > > flatCopy(const In & in,T * out)87 T *flatCopy ( const In &in, T *out ) 88 { 89 for( auto it = in.begin(), end = in.end(); it != end; ++it ) 90 out = flatCopy( *it, out ); 91 return out; 92 } 93 94 95 96 // makeNumPyArray 97 // -------------- 98 99 template< class T, class In > makeNumPyArray(const In & in,const std::vector<std::size_t> & shape)100 pybind11::array_t< T > makeNumPyArray ( const In &in, const std::vector< std::size_t > &shape ) 101 { 102 const std::size_t dim = shape.size(); 103 104 std::size_t size = sizeof( T ); 105 std::vector< std::size_t > stride( dim ); 106 for( std::size_t i = dim; i-- > 0; ) 107 { 108 stride[ i ] = size; 109 size *= shape[ i ]; 110 } 111 112 pybind11::array_t< T > result( pybind11::buffer_info( nullptr, sizeof( T ), pybind11::format_descriptor< T >::value, dim, shape, stride ) ); 113 114 pybind11::buffer_info info = result.request( true ); 115 flatCopy( in, static_cast< T * >( info.ptr ) ); 116 117 return result; 118 } 119 120 121 122 // tessellate 123 // ---------- 124 125 template< class GridView, unsigned int partitions > 126 inline static std::pair< pybind11::array_t< typename GridView::ctype >, pybind11::array_t< int > > tesselate(const GridView & gridView,RefinementIntervals intervals,PartitionSet<partitions> ps)127 tesselate ( const GridView &gridView, RefinementIntervals intervals, PartitionSet< partitions > ps ) 128 { 129 typedef typename GridView::ctype ctype; 130 131 const std::size_t dimGrid = GridView::dimension; 132 const std::size_t dimWorld = GridView::dimensionworld; 133 134 std::vector< FieldVector< ctype, dimWorld > > coords; 135 std::vector< std::array< int, dimGrid+1 > > simplices; 136 137 for( const auto &element : elements( gridView, ps ) ) 138 { 139 const auto &refinement = buildRefinement< dimGrid, double >( element.type(), GeometryTypes::simplex( dimGrid ) ); 140 const std::size_t offset = coords.size(); 141 142 // get coordinates 143 const auto geometry = element.geometry(); 144 for( auto it = refinement.vBegin( intervals ), end = refinement.vEnd( intervals ); it != end; ++it ) 145 coords.push_back( geometry.global( it.coords() ) ); 146 147 // get simplices 148 for( auto it = refinement.eBegin( intervals ), end = refinement.eEnd( intervals ); it != end; ++it ) 149 { 150 std::array< int, dimGrid+1 > simplex; 151 auto indices = it.vertexIndices(); 152 assert( indices.size() == simplex.size() ); 153 std::transform( indices.begin(), indices.end(), simplex.begin(), [ offset ] ( std::size_t i ) { return (i + offset); } ); 154 simplices.push_back( simplex ); 155 } 156 } 157 158 return std::make_pair( makeNumPyArray< ctype >( coords, { coords.size(), dimWorld } ), makeNumPyArray< int >( simplices, { simplices.size(), dimGrid+1 } ) ); 159 } 160 161 template< class GridView, unsigned int partitions > 162 inline static std::pair< pybind11::array_t< typename GridView::ctype >, pybind11::array_t< int > > tesselate(const GridView & gridView,int level,PartitionSet<partitions> ps)163 tesselate ( const GridView &gridView, int level, PartitionSet< partitions > ps ) 164 { 165 return tesselate( gridView, refinementLevels( level ), ps ); 166 } 167 168 template< class GridView, unsigned int partitions > 169 inline static std::vector< pybind11::array_t< typename GridView::ctype > > polygons(const GridView & gridView,PartitionSet<partitions> ps)170 polygons ( const GridView &gridView, PartitionSet< partitions > ps ) 171 { 172 typedef typename GridView::ctype ctype; 173 174 const std::size_t dimWorld = GridView::dimensionworld; 175 176 std::map< std::size_t, std::vector< std::vector< FieldVector< ctype, dimWorld > > > > coords; 177 178 for( const auto &element : elements( gridView, ps ) ) 179 { 180 // get coordinates 181 const auto geometry = element.geometry(); 182 const std::size_t corners = geometry.corners(); 183 std::vector< FieldVector< ctype, dimWorld > > poly( corners ); 184 for( std::size_t i = 0; i != corners; ++i ) 185 poly[ i ] = geometry.corner( i ); 186 187 // Martin: This seems to be limited to two spatial dimensions. 188 if( element.type().isCube() ) 189 std::swap( poly[ 0 ], poly[ 1 ] ); 190 191 coords[ corners ].push_back( poly ); 192 } 193 194 std::vector< pybind11::array_t< typename GridView::ctype > > ret; 195 for( const auto &entry : coords ) 196 ret.push_back( makeNumPyArray< ctype >( entry.second, { entry.second.size(), entry.first, dimWorld } ) ); 197 return ret; 198 } 199 200 template< class GridFunction, unsigned int partitions > 201 inline static 202 // std::pair< 203 // std::vector< pybind11::array_t< typename GridFunction::GridView::ctype > >, 204 // pybind11::array_t< typename GridFunction::GridView::ctype > 205 // > polygonData(const GridFunction & gridFunction,PartitionSet<partitions> ps)206 auto polygonData ( const GridFunction &gridFunction, PartitionSet< partitions > ps ) 207 { 208 typedef typename GridFunction::GridView GridView; 209 typedef typename GridView::ctype ctype; 210 211 const std::size_t dimWorld = GridView::dimensionworld; 212 typedef typename GridFunctionTraits< GridFunction >::Range Range; 213 214 std::map< std::size_t, std::pair< 215 std::vector< std::vector< FieldVector< ctype, dimWorld > > >, 216 std::vector< Range > 217 > > coords; 218 219 const auto &gv = gridView( gridFunction ); 220 auto lf = localFunction( gridFunction ); 221 222 for( const auto &element : elements( gv, ps ) ) 223 { 224 // get coordinates 225 const auto geometry = element.geometry(); 226 const std::size_t corners = geometry.corners(); 227 std::vector< FieldVector< ctype, dimWorld > > poly( corners ); 228 for( std::size_t i = 0; i != corners; ++i ) 229 poly[ i ] = geometry.corner( i ); 230 231 // Martin: This seems to be limited to two spatial dimensions. 232 if( element.type().isCube() ) 233 std::swap( poly[ 0 ], poly[ 1 ] ); 234 235 lf.bind(element); 236 coords[ corners ].first.push_back( poly ); 237 // for polygons we can't use a reference element but one could use 238 // the following for elements with type not none: 239 // auto ref = ReferenceElements< typename GridView::ctype, dimGrid >::general( element.type() ); 240 // coords[ corners ].second.push_back( lf( ref.position(0,0) ) ); 241 coords[ corners ].second.push_back( lf( geometry.local( geometry.center() ) ) ); 242 lf.unbind(); 243 } 244 245 std::vector< pybind11::array_t< typename GridView::ctype > > ret; 246 std::vector< pybind11::array_t< typename GridView::ctype > > values; 247 for( const auto &entry : coords ) 248 { 249 ret.push_back( makeNumPyArray< ctype >( entry.second.first, { entry.second.first.size(), entry.first, dimWorld } ) ); 250 values.push_back( Python::makeNumPyArray< typename FieldTraits< Range >::field_type >( entry.second.second, { entry.second.second.size(), GetDimension<Range>::value } ) ); 251 } 252 return std::make_pair(ret, values); 253 } 254 255 template< class GridView, unsigned int partitions > 256 inline static std::pair< pybind11::array_t< typename GridView::ctype >, pybind11::array_t< int > > tesselate(const GridView & gridView,PartitionSet<partitions> ps)257 tesselate ( const GridView &gridView, PartitionSet< partitions > ps ) 258 { 259 return tesselate( gridView, 0, ps ); 260 } 261 262 template< class GridView > 263 inline static std::pair< pybind11::array_t< typename GridView::ctype >, pybind11::array_t< int > > tesselate(const GridView & gridView,int level=0)264 tesselate ( const GridView &gridView, int level = 0 ) 265 { 266 return tesselate( gridView, level, Partitions::all ); 267 } 268 269 template< class GridView > 270 inline static std::vector< pybind11::array_t< typename GridView::ctype > > polygons(const GridView & gridView)271 polygons ( const GridView &gridView ) 272 { 273 return polygons( gridView, Partitions::all ); 274 } 275 276 template< class GridFunction > 277 inline static 278 // std::pair< 279 // std::vector< pybind11::array_t< typename GridFunction::GridView::ctype > >, 280 // pybind11::array_t< typename GridFunction::GridView::ctype > 281 // > polygonData(const GridFunction & gridFunction)282 auto polygonData ( const GridFunction &gridFunction ) 283 { 284 return polygonData( gridFunction, Partitions::all ); 285 } 286 287 // pointData 288 // --------- 289 290 template< class GridFunction, unsigned int partitions > pointData(const GridFunction & gridFunction,int level,PartitionSet<partitions> ps)291 inline static auto pointData ( const GridFunction &gridFunction, int level, PartitionSet< partitions > ps ) 292 { 293 typedef typename GridFunctionTraits< GridFunction >::Element Element; 294 typedef typename GridFunctionTraits< GridFunction >::LocalCoordinate LocalCoordinate; 295 typedef typename GridFunctionTraits< GridFunction >::Range Range; 296 297 typedef typename FieldTraits< LocalCoordinate >::field_type ctype; 298 299 const auto &gv = gridView( gridFunction ); 300 301 auto lf = localFunction( gridFunction ); 302 std::vector< Range > values; 303 auto refLevel = refinementLevels(level); 304 for( const auto &element : elements( gv, ps ) ) 305 { 306 lf.bind( element ); 307 const auto &refinement = buildRefinement< Element::mydimension, ctype >( element.type(), GeometryTypes::simplex( Element::dimension ) ); 308 for( auto it = refinement.vBegin( refLevel ), end = refinement.vEnd( refLevel ); it != end; ++it ) 309 values.push_back( lf( it.coords() ) ); 310 lf.unbind(); 311 } 312 313 return Python::makeNumPyArray< typename FieldTraits< Range >::field_type >( values, { values.size(), GetDimension<Range>::value } ); 314 } 315 316 template< class GridFunction, unsigned int partitions > pointData(const GridFunction & gridFunction,PartitionSet<partitions> ps)317 inline static auto pointData ( const GridFunction &gridFunction, PartitionSet< partitions > ps ) 318 { 319 return pointData( gridFunction, 0, ps ); 320 } 321 322 template< class GridFunction > pointData(const GridFunction & gridFunction,int level=0)323 inline static auto pointData ( const GridFunction &gridFunction, int level = 0 ) 324 { 325 return pointData( gridFunction, level, Partitions::all ); 326 } 327 328 329 330 // cellData 331 // -------- 332 333 template< class GridFunction, unsigned int partitions > cellData(const GridFunction & gridFunction,int level,PartitionSet<partitions> ps)334 inline static auto cellData ( const GridFunction &gridFunction, int level, PartitionSet< partitions > ps ) 335 { 336 typedef typename GridFunctionTraits< GridFunction >::Element Element; 337 typedef typename GridFunctionTraits< GridFunction >::LocalCoordinate LocalCoordinate; 338 typedef typename GridFunctionTraits< GridFunction >::Range Range; 339 340 typedef typename FieldTraits< LocalCoordinate >::field_type ctype; 341 342 const auto &gv = gridView( gridFunction ); 343 344 auto lf = localFunction( gridFunction ); 345 std::vector< Range > values; 346 auto refLevel = refinementLevels(0); 347 for( const Element &element : entities( gv, Dune::Codim< Element::codimension >(), ps ) ) 348 { 349 lf.bind( element ); 350 const auto &refinement = buildRefinement< Element::mydimension, ctype >( element.type(), GeometryTypes::simplex( Element::mydimension ) ); 351 for( auto it = refinement.eBegin( refLevel ), end = refinement.eEnd( refLevel ); it != end; ++it ) 352 values.push_back( lf( it.coords() ) ); 353 lf.unbind(); 354 } 355 356 return Python::makeNumPyArray< typename FieldTraits< Range >::field_type >( values, { values.size(), GetDimension<Range>::value } ); 357 } 358 359 template< class GridFunction, unsigned int partitions > cellData(const GridFunction & gridFunction,PartitionSet<partitions> ps)360 inline static auto cellData ( const GridFunction &gridFunction, PartitionSet< partitions > ps ) 361 { 362 return cellData( gridFunction, 0, ps ); 363 } 364 365 template< class GridFunction > cellData(const GridFunction & gridFunction,int level=0)366 inline static auto cellData ( const GridFunction &gridFunction, int level = 0 ) 367 { 368 return cellData( gridFunction, level, Partitions::all ); 369 } 370 371 } // namespace FemPy 372 373 } // namespace Dune 374 375 #endif // #ifndef DUNE_PYTHON_GRID_NUMPY_HH 376