1#ifndef DUNE_FEM_SPACE_FOURIER_INTERPOLATE_HH
2#define DUNE_FEM_SPACE_FOURIER_INTERPOLATE_HH
7#include <dune/common/exceptions.hh>
8#include <dune/common/fmatrix.hh>
9#include <dune/common/fvector.hh>
11#include <dune/grid/common/partitionset.hh>
12#include <dune/grid/common/rangegenerators.hh>
28 template<
class Gr
idFunction,
class DiscreteFunction,
unsigned int partitions >
29 static inline std::enable_if_t< std::is_convertible< GridFunction, HasLocalFunction >::value && IsFourierDiscreteFunctionSpace< typename DiscreteFunction::DiscreteFunctionSpaceType >::value >
30 interpolate (
const GridFunction &u, DiscreteFunction &v, PartitionSet< partitions > ps )
32 typedef typename DiscreteFunction::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
34 typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
35 typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
36 typedef typename DiscreteFunctionSpaceType::RangeFieldType RangeFieldType;
38 const int dimDomain = DiscreteFunctionSpaceType::dimDomain;
39 const int dimRange = DiscreteFunctionSpaceType::dimRange;
44 FieldMatrix< RangeFieldType, size, size > mat( 0 );
45 FieldMatrix< RangeFieldType, dimRange, size > rhs( 0 );
47 const auto &functionSet = v.space().functionSet();
48 if( functionSet.order() == std::numeric_limits< int >::max() )
49 DUNE_THROW( InvalidStateException,
"Cannot interpolate to Fourier space if quadrature order is not specified" );
50 typename GridFunction::LocalFunctionType uLocal( u );
56 for(
const auto entity : elements( v.gridPart(), Partitions::interior ) )
58 const auto geometry = entity.geometry();
60 uLocal.init( entity );
65 const auto x = geometry.global( qp.position() );
66 const RangeFieldType weight = geometry.integrationElement( x ) * qp.weight();
69 FieldVector< RangeFieldType, size > values;
70 functionSet.evaluateEach( x, [ &values ] ( std::size_t i,
const RangeFieldType &v ) { values[ i ] = v; } );
73 for(
int i = 0; i < size; ++i )
74 mat[ i ].
axpy( values[ i ]*weight, values );
78 uLocal.evaluate( qp, uValue );
81 for(
int i = 0; i < dimRange; ++i )
82 rhs[ i ].
axpy( uValue[ i ]*weight, values );
87 std::array< RangeFieldType, (size+dimRange)*size > buffer;
88 for(
int i = 0; i < size; ++i )
89 std::copy_n( mat[ i ].begin(), size, buffer.begin() + i*size );
90 for(
int i = 0; i < dimRange; ++i )
91 std::copy_n( rhs[ i ].begin(), size, buffer.begin() + (i+size)*size );
92 v.gridPart().comm().sum( buffer.data(), buffer.size() );
93 for(
int i = 0; i < size; ++i )
94 std::copy_n( buffer.begin() + i*size, size, mat[ i ].begin() );
95 for(
int i = 0; i < dimRange; ++i )
96 std::copy_n( buffer.begin() + (i+size)*size, size, rhs[ i ].begin() );
100 FieldMatrix< RangeFieldType, dimRange, size >
dofs( 0 );
101 for(
int i = 0; i < dimRange; ++i )
102 mat.umv( rhs[ i ],
dofs[ i ] );
105 for(
int i = 0; i < size; ++i )
106 for(
int j = 0; j < dimRange; ++j )
107 v.dofVector()[ 0 ][ i*dimRange + j ] =
dofs[ j ][ i ];
static void interpolate(const GridFunction &u, DiscreteFunction &v)
perform native interpolation of a discrete function space
Definition: common/interpolate.hh:54
Definition: bindguard.hh:11
IteratorRange< typename DF::DofIteratorType > dofs(DF &df)
Iterates over all DOFs.
Definition: rangegenerators.hh:76
void axpy(const T &a, const T &x, T &y)
Definition: space/basisfunctionset/functor.hh:38
Definition: space/fourier/functionset.hh:28