dune-fem  2.8-git
discretefunction_inline.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_DISCRETEFUNCTION_INLINE_HH
2 #define DUNE_FEM_DISCRETEFUNCTION_INLINE_HH
3 
4 #include <fstream>
5 
6 #include <dune/geometry/referenceelements.hh>
7 
13 
14 #include "discretefunction.hh"
15 
16 namespace Dune
17 {
18 
19  namespace Fem
20  {
21 
22  // DiscreteFunctionDefault
23  // -----------------------
24 
25  template< class Impl >
27  :: DiscreteFunctionDefault ( const std::string &name,
28  const DiscreteFunctionSpaceType &dfSpace )
29  : dfSpace_( referenceToSharedPtr( dfSpace ) ),
30  ldvStack_( std::max( std::max( sizeof( DofType ), sizeof( DofType* ) ),
31  sizeof(typename LocalDofVectorType::value_type) ) // for PetscDiscreteFunction
32  * space().blockMapper().maxNumDofs() * DiscreteFunctionSpaceType::localBlockSize ),
33  ldvAllocator_( &ldvStack_ ),
34  localFunction_( space() ),
35  name_( name ),
36  scalarProduct_( dfSpace )
37  {
38  }
39 
40 
41  template< class Impl >
42  inline DiscreteFunctionDefault< Impl >::DiscreteFunctionDefault ( std::string name, std::shared_ptr< const DiscreteFunctionSpaceType > dfSpace )
43  : dfSpace_( std::move( dfSpace ) ),
44  ldvStack_( std::max( std::max( sizeof( DofType ), sizeof( DofType* ) ),
45  sizeof(typename LocalDofVectorType::value_type) ) // for PetscDiscreteFunction
46  * space().blockMapper().maxNumDofs() * DiscreteFunctionSpaceType::localBlockSize ),
47  ldvAllocator_( &ldvStack_ ),
48  localFunction_( space() ),
49  name_( std::move( name ) ),
50  scalarProduct_( dfSpace )
51  {}
52 
53 
54  template< class Impl >
56  : BaseType( static_cast< const BaseType & >( other ) ),
57  dfSpace_( other.dfSpace_ ),
58  ldvStack_( std::max( std::max( sizeof( DofType ), sizeof( DofType* ) ),
59  sizeof(typename LocalDofVectorType::value_type) ) // for PetscDiscreteFunction
60  * space().blockMapper().maxNumDofs() * DiscreteFunctionSpaceType::localBlockSize ),
61  ldvAllocator_( &ldvStack_ ),
62  localFunction_( space() ),
63  name_( other.name_ ),
64  scalarProduct_( other.scalarProduct_ )
65  {
66  if( other.assembleOperation_ != std::type_index( typeid( void ) ) )
67  DUNE_THROW( InvalidStateException, "Cannot copy discrete function during assembly" );
68  assert( other.assembleCount_ == 0 );
69  }
70 
71 
72  template< class Impl >
74  :: DiscreteFunctionDefault ( DiscreteFunctionDefault && other )
75  : BaseType( static_cast< BaseType&& >( other ) ),
76  dfSpace_( std::move( other.dfSpace_ ) ),
77  ldvStack_( std::move( other.ldvStack_ ) ),
78  ldvAllocator_( &ldvStack_ ),
79  localFunction_( space() ),
80  name_( std::move( other.name_ ) ),
81  scalarProduct_( std::move( other.scalarProduct_ ) )
82  {
83  if( other.assembleOperation_ != std::type_index( typeid( void ) ) )
84  DUNE_THROW( InvalidStateException, "Cannot move discrete function during assembly" );
85  assert( other.assembleCount_ == 0 );
86  }
87 
88 
89  template< class Impl >
91  :: print ( std::ostream &out ) const
92  {
93  const auto end = BaseType :: dend();
94  for( auto dit = BaseType :: dbegin(); dit != end; ++dit )
95  out << (*dit) << std::endl;
96  }
97 
98 
99  template< class Impl >
102  {
103  const auto end = BaseType :: dend();
104  for( auto it = BaseType :: dbegin(); it != end; ++it )
105  {
106  if( ! std::isfinite( *it ) )
107  return false ;
108  }
109 
110  return true;
111  }
112 
113 
114  template< class Impl >
115  template< class DFType >
118  {
119  if( BaseType::size() != g.size() )
120  DUNE_THROW(InvalidStateException,"DiscreteFunctionDefault: sizes do not match in axpy");
121 
122  // apply axpy to all dofs from g
123  const auto end = BaseType::dend();
124  auto git = g.dbegin();
125  for( auto it = BaseType::dbegin(); it != end; ++it, ++git )
126  *it += s * (*git );
127  }
128 
129 
130  template< class Impl >
131  template< class DFType >
134  {
135  if( BaseType::size() != g.size() )
136  {
137  std::cout << BaseType::size() << " vs " << g.size() << std::endl;
138  DUNE_THROW(InvalidStateException,"DiscreteFunctionDefault: sizes do not match in assign");
139  }
140 
141  // copy all dofs from g to this
142  const auto end = BaseType::dend();
143  auto git = g.dbegin();
144  for( auto it = BaseType::dbegin(); it != end; ++it, ++git )
145  *it = *git;
146  }
147 
148 
149  template< class Impl >
150  template< class Operation >
151  inline typename DiscreteFunctionDefault< Impl >
152  :: template CommDataHandle< Operation > :: Type
153  DiscreteFunctionDefault< Impl > :: dataHandle ( const Operation &operation )
154  {
155  return BaseType :: space().createDataHandle( asImp(), operation );
156  }
157 
158 
159  template< class Impl >
160  template< class Functor >
162  ::evaluateGlobal ( const DomainType &x, Functor functor ) const
163  {
164  typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
165  EntitySearch< GridPartType, EntityType::codimension > entitySearch( BaseType::space().gridPart() );
166 
167  const auto& entity(entitySearch( x ));
168  const auto geometry = entity.geometry();
169 
170  // bind local function to entity
171  auto guard = bindGuard( localFunction_, entity );
172  // obtain dofs (since it's a temporary local function)
173  getLocalDofs( entity, localFunction_.localDofVector());
174  // evaluate functor
175  functor( geometry.local( x ), localFunction_ );
176  }
177 
178 
179  template< class Impl >
180  template< class DFType >
181  inline typename DiscreteFunctionDefault< Impl > :: DiscreteFunctionType &
184  {
185  if( BaseType::size() != g.size() )
186  DUNE_THROW(InvalidStateException,"DiscreteFunctionDefault: sizes do not match in operator +=");
187 
188  const auto end = BaseType::dend();
189  auto git = g.dbegin();
190  for( auto it = BaseType::dbegin(); it != end; ++it, ++git )
191  *it += *git;
192  return asImp();
193  }
194 
195 
196  template< class Impl >
197  template< class DFType >
198  inline typename DiscreteFunctionDefault< Impl > :: DiscreteFunctionType &
201  {
202  if( BaseType::size() != g.size() )
203  DUNE_THROW(InvalidStateException,"DiscreteFunctionDefault: sizes do not match in operator -=");
204 
205  const auto end = BaseType :: dend();
206  auto git = g.dbegin();
207  for( auto it = BaseType :: dbegin(); it != end; ++it, ++git )
208  *it -= *git;
209  return asImp();
210  }
211 
212 
213  template< class Impl >
214  template< class StreamTraits >
217  {
218  auto versionId = in.readUnsignedInt();
219  if( versionId < DUNE_VERSION_ID(0,9,1) )
220  DUNE_THROW( IOError, "Trying to read outdated file." );
221  else if( versionId > DUNE_MODULE_VERSION_ID(DUNE_FEM) )
222  std :: cerr << "Warning: Reading discrete function from newer version: "
223  << versionId << std :: endl;
224 
225  // verify space id for files written with dune-fem version 1.5 or newer
226  if( versionId >= DUNE_VERSION_ID(1,5,0) )
227  {
228  // make sure that space of discrete function matches the space
229  // of the data that was written
230  const auto spaceId = space().type();
231  int mySpaceIdInt;
232  in >> mySpaceIdInt;
233  const auto mySpaceId = static_cast<DFSpaceIdentifier>(mySpaceIdInt);
234 
235  if( spaceId != mySpaceId )
236  DUNE_THROW( IOError, "Trying to read discrete function from different space: DFSpace (" << spaceName( spaceId ) << ") != DataSpace (" << spaceName( mySpaceId ) << ")" );
237  }
238 
239  // read name
240  in >> name_;
241 
242  // read size as integer
243  int32_t mysize;
244  in >> mysize;
245 
246  // check size
247  if( static_cast< int >( mysize ) != BaseType::size() )
248  DUNE_THROW( IOError, "Trying to read discrete function of different size." );
249  if( BaseType::size() != static_cast< int >( this->space().size() ) )
250  DUNE_THROW( InvalidStateException, "Trying to read discrete function in uncompressed state." );
251 
252  // read all dofs
253  typedef typename DiscreteFunctionSpaceType::BlockMapperType BlockMapperType;
254  typedef typename DiscreteFunctionSpaceType::LocalBlockIndices LocalBlockIndices;
255  typedef typename BlockMapperType::SizeType SizeType;
256 
257  const BlockMapperType &blockMapper = this->space().blockMapper();
258  const SizeType nBlocks = blockMapper.size();
259  for( SizeType i = 0; i < nBlocks; ++i )
260  {
261  auto &&block = this->dofVector()[ i ];
262  Hybrid::forEach( LocalBlockIndices(), [ &in, &block ] ( auto j ) { in >> block[ j ]; } );
263  }
264 
265  // a discrete function that is part of backup/restore
266  // most likely needs to keep the dofs consistent
267  asImp().enableDofCompression();
268  }
269 
270 
271  template< class Impl >
272  template< class StreamTraits >
275  {
276  unsigned int versionId = DUNE_MODULE_VERSION_ID(DUNE_FEM);
277  out << versionId ;
278 
279  // write space id to for testing when function is read
280  auto spaceId = space().type();
281  out << spaceId ;
282 
283  // write name
284  out << name_;
285 
286  // only allow write when vector is compressed
287  if( BaseType::size() != static_cast< int >( this->space().size() ) )
288  DUNE_THROW( InvalidStateException, "Trying to write DiscreteFunction in uncompressed state." );
289 
290  // write size as 32-bit integer
291  const int32_t mysize = BaseType::size();
292  out << mysize;
293 
294  // write all dofs
295  typedef typename DiscreteFunctionSpaceType::BlockMapperType BlockMapperType;
296  typedef typename DiscreteFunctionSpaceType::LocalBlockIndices LocalBlockIndices;
297  typedef typename BlockMapperType::SizeType SizeType;
298 
299  const BlockMapperType &blockMapper = this->space().blockMapper();
300  const SizeType nBlocks = blockMapper.size();
301  for( SizeType i = 0; i < nBlocks; ++i )
302  {
303  const auto block = this->dofVector()[ i ];
304  Hybrid::forEach( LocalBlockIndices(), [ &out, &block ] ( auto j ) { out << block[ j ]; } );
305  }
306  }
307 
308 
309  template< class Impl >
312  {
313  typedef typename DiscreteFunctionSpaceType::IndexSetType IndexSetType;
314  IndexSetType& indexSet = (IndexSetType&)space().indexSet();
316  {
317  auto persistentIndexSet = Dune::Fem::Capabilities::isPersistentIndexSet< IndexSetType >::map( indexSet );
318 
319  // this marks the index set in the DofManager's list of index set as persistent
320  if( persistentIndexSet )
321  persistentIndexSet->addBackupRestore();
322  }
323  }
324 
325  template< class Impl >
328  {
329  typedef typename DiscreteFunctionSpaceType::IndexSetType IndexSetType;
330  IndexSetType& indexSet = (IndexSetType&)space().indexSet();
332  {
333  auto persistentIndexSet = Dune::Fem::Capabilities::isPersistentIndexSet< IndexSetType >::map( indexSet );
334 
335  // this unmarks the index set in the DofManager's list of index set as persistent
336  if( persistentIndexSet )
337  persistentIndexSet->removeBackupRestore();
338  }
339  }
340 
341 
342  template< class Impl >
343  template< class DFType >
346  {
347  if( BaseType :: size() != g.size() )
348  return false;
349 
350  const auto end = BaseType :: dend();
351 
352  auto fit = BaseType :: dbegin();
353  auto git = g.dbegin();
354  for( ; fit != end; ++fit, ++git )
355  {
356  if( std::abs( *fit - *git ) > 1e-15 )
357  return false;
358  }
359 
360  return true;
361  }
362 
363 
364  // Stream Operators
365  // ----------------
366 
375  template< class Impl >
376  inline std :: ostream &
377  operator<< ( std :: ostream &out,
379  {
380  df.print( out );
381  return out;
382  }
383 
384 
385 
395  template< class StreamTraits, class Impl >
399  {
400  df.write( out );
401  return out;
402  }
403 
404 
405 
415  template< class StreamTraits, class Impl >
416  inline InStreamInterface< StreamTraits > &
419  {
420  df.read( in );
421  return in;
422  }
423 
424  } // end namespace Fem
425 
426 } // end namespace Dune
427 #endif // #ifndef DUNE_FEM_DISCRETEFUNCTION_INLINE_HH
DFSpaceIdentifier
enumerator for identification of spaces
Definition: discretefunctionspace.hh:94
std::string spaceName(const DFSpaceIdentifier id)
Definition: discretefunctionspace.hh:109
Definition: bindguard.hh:11
OutStreamInterface< StreamTraits > & operator<<(OutStreamInterface< StreamTraits > &out, const DiscreteFunctionInterface< Impl > &df)
write a discrete function into an output stream
Definition: discretefunction_inline.hh:397
Double abs(const Double &a)
Definition: double.hh:871
static double max(const Double &v, const double p)
Definition: double.hh:398
static auto bindGuard(Object &object, Args &&... args) -> std::enable_if_t< isBindable< Object, Args... >::value, BindGuard< Object > >
Definition: bindguard.hh:67
static std::shared_ptr< T > referenceToSharedPtr(T &t)
Definition: memory.hh:19
InStreamInterface< StreamTraits > & operator>>(InStreamInterface< StreamTraits > &in, DiscreteFunctionInterface< Impl > &df)
read a discrete function from an input stream
Definition: discretefunction_inline.hh:417
static void forEach(IndexRange< T, sz > range, F &&f)
Definition: hybrid.hh:129
Definition: common/discretefunction.hh:578
void print(std ::ostream &out) const
print all DoFs to a stream (for debugging purposes)
Definition: discretefunction_inline.hh:91
DiscreteFunctionDefault(const std::string &name, const DiscreteFunctionSpaceType &dfSpace)
Constructor storing discrete function space and local function factory.
Definition: discretefunction_inline.hh:27
DiscreteFunctionType & operator-=(const DiscreteFunctionInterface< DFType > &g)
substract all degrees of freedom from given discrete function using the dof iterators
Definition: discretefunction_inline.hh:200
void write(OutStreamInterface< StreamTraits > &out) const
write the discrete function into a stream
Definition: discretefunction_inline.hh:274
DiscreteFunctionType & operator+=(const DiscreteFunctionInterface< DFType > &g)
add another discrete function to this one
Definition: discretefunction_inline.hh:183
bool dofsValid() const
check for NaNs
Definition: discretefunction_inline.hh:101
DofVectorType::SizeType SizeType
size type of the block vector
Definition: common/discretefunction.hh:646
std::type_index assembleOperation_
Definition: common/discretefunction.hh:1042
void evaluateGlobal(const DomainType &x, Functor functor) const
evaluate functor in global coordinate
Definition: discretefunction_inline.hh:162
void read(InStreamInterface< StreamTraits > &in)
read the discrete function from a stream
Definition: discretefunction_inline.hh:216
virtual void insertSubData()
Definition: discretefunction_inline.hh:311
BaseType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType
type of discrete function space
Definition: common/discretefunction.hh:600
void assign(const DiscreteFunctionInterface< DFType > &g)
Definition: discretefunction_inline.hh:133
virtual void removeSubData()
Definition: discretefunction_inline.hh:327
CommDataHandle< Operation >::Type dataHandle(const Operation &operation)
return reference to data handle object
Traits ::LocalDofVectorType LocalDofVectorType
type of LocalDofVector
Definition: common/discretefunction.hh:628
void axpy(const RangeFieldType &s, const DiscreteFunctionInterface< DFType > &g)
axpy operation
Definition: discretefunction_inline.hh:117
std::size_t assembleCount_
Definition: common/discretefunction.hh:1043
bool operator==(const DiscreteFunctionInterface< DFType > &g) const
Definition: discretefunction_inline.hh:345
Definition: common/discretefunction.hh:86
DiscreteFunctionSpaceType ::RangeFieldType RangeFieldType
type of range field, i.e. dof type
Definition: common/discretefunction.hh:109
void write(OutStreamInterface< StreamTraits > &out) const
write the discrete function into a stream
Definition: common/discretefunction.hh:537
ConstDofIteratorType dbegin() const
obtain an iterator pointing to the first DoF (read-only)
Definition: common/discretefunction.hh:354
void print(std ::ostream &out) const
print all DoFs to a stream (for debugging purposes)
Definition: common/discretefunction.hh:437
void read(InStreamInterface< StreamTraits > &in)
read the discrete function from a stream
Definition: common/discretefunction.hh:527
Traits ::DofType DofType
Definition: common/discretefunction.hh:136
DiscreteFunctionSpaceType ::DomainType DomainType
type of domain, i.e. type of coordinates
Definition: common/discretefunction.hh:111
DiscreteFunctionSpaceType::GridPartType GridPartType
type of the underlying grid part
Definition: common/discretefunction.hh:118
int size() const
obtain total number of DoFs
Definition: common/discretefunction.hh:333
Definition: entitysearch.hh:131
capability for persistent index sets
Definition: persistentindexset.hh:92
static constexpr PersistentIndexSetInterface * map(IndexSet &indexSet) noexcept
please doc me
Definition: persistentindexset.hh:101
abstract interface for an output stream
Definition: streams.hh:46
abstract interface for an input stream
Definition: streams.hh:179