File: bvector.hh

package info (click to toggle)
dune-istl 2.11.0-1
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 3,504 kB
  • sloc: cpp: 34,844; python: 182; sh: 3; makefile: 3
file content (236 lines) | stat: -rw-r--r-- 9,712 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_PYTHON_ISTL_BVECTOR_HH
#define DUNE_PYTHON_ISTL_BVECTOR_HH

#include <cstddef>

#include <stdexcept>
#include <string>
#include <type_traits>
#include <utility>

#include <dune/common/typeutilities.hh>

//this added otherwise insert class wasn't possible on line ~190
#include <dune/python/common/typeregistry.hh>
#include <dune/python/common/fvecmatregistry.hh>
#include <dune/python/common/string.hh>
#include <dune/python/common/vector.hh>
#include <dune/python/istl/iterator.hh>
#include <dune/python/pybind11/operators.h>
#include <dune/python/pybind11/pybind11.h>

#include <dune/istl/bvector.hh>
#include <dune/istl/blocklevel.hh>

namespace Dune
{

  namespace Python
  {

    namespace detail
    {

      template< class K, int n >
      inline static void copy ( const char *ptr, const ssize_t *shape, const ssize_t *strides, Dune::FieldVector< K, n > &v )
      {
        if( *shape != static_cast< ssize_t >( n ) )
          throw pybind11::value_error( "Invalid buffer size: " + std::to_string( *shape ) + " (should be: " + std::to_string( n ) + ")." );

        for( ssize_t i = 0; i < static_cast< ssize_t >( n ); ++i )
          v[ i ] = *reinterpret_cast< const K * >( ptr + i*(*strides) );
      }


      template< class B, class A >
      inline static void copy ( const char *ptr, const ssize_t *shape, const ssize_t *strides, Dune::BlockVector< B, A > &v )
      {
        v.resize( *shape );
        for( ssize_t i = 0; i < *shape; ++i )
          copy( ptr + i*(*strides), shape+1, strides+1, v[ i ] );
      }


      template< class BlockVector >
      inline static void copy ( pybind11::buffer buffer, BlockVector &v )
      {
        typedef typename BlockVector::field_type field_type;
        typedef typename BlockVector::size_type  size_type;

        pybind11::buffer_info info = buffer.request();

        if( info.format != pybind11::format_descriptor< field_type >::format() )
          throw pybind11::value_error( "Incompatible buffer format." );
        if( size_type(info.ndim) != blockLevel<BlockVector>() )
          throw pybind11::value_error( "Block vectors can only be initialized from one-dimensional buffers." );

        copy( static_cast< const char * >( info.ptr ), info.shape.data(), info.strides.data(), v );
      }



      // blockVectorGetItem
      // ------------------

      template< class BlockVector >
      inline static pybind11::object blockVectorGetItem ( const pybind11::object &vObj, BlockVector &v, typename BlockVector::size_type index )
      {
        auto pos = v.find( index );
        if( pos == v.end() )
          throw pybind11::index_error( "Index " + std::to_string( index ) + " does not exist in block vector." );
        pybind11::object result = pybind11::cast( *pos, pybind11::return_value_policy::reference );
        pybind11::detail::keep_alive_impl( result, vObj );
        return result;
      }

    } // namespace detail



    // to_string
    // ---------

    template< class X >
    inline static auto to_string ( const X &x )
      -> std::enable_if_t< std::is_base_of< Imp::block_vector_unmanaged< typename X::block_type, typename X::size_type >, X >::value, std::string >
    {
      return "(" + join( ", ", [] ( auto &&x ) { return to_string( x ); }, x.begin(), x.end() ) + ")";
    }



    // registserBlockVector
    // --------------------

    template< class BlockVector, class... options >
    inline void registerBlockVector ( pybind11::class_< BlockVector, options... > cls )
    {
      typedef typename BlockVector::field_type field_type;
      typedef typename BlockVector::block_type block_type;
      typedef typename BlockVector::size_type size_type;

      using pybind11::operator""_a;

      cls.def( "assign", [] ( BlockVector &self, const BlockVector &x ) { self = x; }, "x"_a );

      cls.def( "copy", [] ( const BlockVector &self ) { return new BlockVector( self ); } );

      cls.def( "__getitem__", [] ( const pybind11::object &self, size_type index ) {
          return detail::blockVectorGetItem( self, pybind11::cast< BlockVector & >( self ), index );
        } );
      cls.def( "__getitem__", [] ( const pybind11::object &self, pybind11::iterable index ) {
          BlockVector &v = pybind11::cast< BlockVector & >( self );
          pybind11::tuple refs( pybind11::len( index ) );
          std::size_t j = 0;
          for( pybind11::handle i : index )
            refs[ j++ ] = detail::blockVectorGetItem( self, v, pybind11::cast< size_type >( i ) );
          return refs;
        } );

      cls.def( "__setitem__", [] ( BlockVector &self, size_type index, block_type value ) {
          auto pos = self.find( index );
          if( pos != self.end() )
            *pos = value;
          else
            throw pybind11::index_error();
        } );
      cls.def( "__setitem__", [] ( BlockVector &self, pybind11::slice index, pybind11::iterable value ) {
          std::size_t start, stop, step, length;
          index.compute( self.N(), &start, &stop, &step, &length );
          for( auto v : value )
          {
            if( start >= stop )
              throw pybind11::value_error( "too many values passed" );
            auto pos = self.find( start );
            if( pos != self.end() )
              *pos = pybind11::cast< block_type >( v );
            else
              throw pybind11::index_error();
            start += step;
          }
          if( start < stop )
            throw pybind11::value_error( "too few values passed" );
        } );

      cls.def( "__len__", [] ( const BlockVector &self ) { return self.N(); } );

      detail::registerOneTensorInterface( cls );
      detail::registerISTLIterators( cls );

      cls.def( "__iadd__", [] ( BlockVector &self, const BlockVector& x ) -> BlockVector & { self += x; return self; } );
      cls.def( "__isub__", [] ( BlockVector &self, const BlockVector& x ) -> BlockVector & { self -= x; return self; } );
      cls.def( "__imul__", [] ( BlockVector &self, field_type x ) -> BlockVector & { self *= x; return self; } );
      cls.def( "__idiv__", [] ( BlockVector &self, field_type x ) -> BlockVector & { self /= x; return self; } );
      cls.def( "__itruediv__", [] ( BlockVector &self, field_type x ) -> BlockVector & { self /= x; return self; } );

      cls.def( "__add__", [] ( const BlockVector &self, const BlockVector &x ) { BlockVector *copy = new BlockVector( self ); *copy += x; return copy; } );
      cls.def( "__sub__", [] ( const BlockVector &self, const BlockVector &x ) { BlockVector *copy = new BlockVector( self ); *copy -= x; return copy; } );

      cls.def( "__div__", [] ( const BlockVector &self, field_type x ) { BlockVector *copy = new BlockVector( self ); *copy /= x; return copy; } );
      cls.def( "__truediv__", [] ( const BlockVector &self, field_type x ) { BlockVector *copy = new BlockVector( self ); *copy /= x; return copy; } );
      cls.def( "__mul__", [] ( const BlockVector &self, field_type x ) { BlockVector *copy = new BlockVector( self ); *copy *= x; return copy; } );
      cls.def( "__rmul__", [] ( const BlockVector &self, field_type x ) { BlockVector *copy = new BlockVector( self ); *copy *= x; return copy; } );
    }



    // registerBlockVector
    // --------------------

    //for the new bindings and arbitrary block size haven't
    //the generator actually takes the scope into account which is why we do nothing with it here
    //so when doing a dune.istl blockvector it doesn't actually define any of the rest of the bindings
    template< class BlockVector, class ... options >
    void registerBlockVector ( pybind11::handle /*scope*/, pybind11::class_<BlockVector, options ... > cls )
    {
      typedef typename BlockVector::size_type size_type;
      using pybind11::operator""_a;

      registerBlockVector( cls );

      cls.def( pybind11::init( [] () { return new BlockVector(); } ) );
      cls.def( pybind11::init( [] ( size_type size ) { return new BlockVector( size ); } ), "size"_a );

      cls.def( pybind11::init( [] ( pybind11::buffer buffer ) {
          BlockVector *self = new BlockVector();
          detail::copy( buffer, *self );
          return self;
        } ) );


      // cls.def( "__str__", [] ( const BlockVector &self ) { return to_string( self ); } );

      cls.def( "assign", [] ( BlockVector &self, pybind11::buffer buffer ) { detail::copy( buffer, self ); }, "buffer"_a );

      cls.def_property_readonly( "capacity", [] ( const BlockVector &self ) { return self.capacity(); } );

      cls.def( "resize", [] ( BlockVector &self, size_type size ) { self.resize( size ); }, "size"_a );

    }

    //the auto class is needed so that run.algorithm can properly work
    template< class BlockVector >
    inline pybind11::class_< BlockVector > registerBlockVector ( pybind11::handle scope, const char *clsName = "BlockVector" )
    {
      //typedef typename BlockVector::size_type size_type;

      using pybind11::operator""_a;

      int rows = BlockVector::block_type::dimension;
      std::string vectorTypename = "Dune::BlockVector< Dune::FieldVector< double, "+ std::to_string(rows) + " > >";
      auto cls = Dune::Python::insertClass< BlockVector >( scope, clsName, Dune::Python::GenerateTypeName(vectorTypename), Dune::Python::IncludeFiles{"dune/istl/bvector.hh","dune/python/istl/bvector.hh"});

      if (cls.second)
        registerBlockVector( scope, cls.first );
      return cls.first;
    }



  } // namespace Python

} // namespace Dune

#endif // #ifndef DUNE_PYTHON_ISTL_BVECTOR_HH