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