dune-fem  2.8-git
fourier/interpolate.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_SPACE_FOURIER_INTERPOLATE_HH
2 #define DUNE_FEM_SPACE_FOURIER_INTERPOLATE_HH
3 
4 #include <limits>
5 #include <type_traits>
6 
7 #include <dune/common/exceptions.hh>
8 #include <dune/common/fmatrix.hh>
9 #include <dune/common/fvector.hh>
10 
11 #include <dune/grid/common/partitionset.hh>
12 #include <dune/grid/common/rangegenerators.hh>
13 
18 
19 namespace Dune
20 {
21 
22  namespace Fem
23  {
24 
25  // interpolate
26  // -----------
27 
28  template< class GridFunction, 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 )
31  {
32  typedef typename DiscreteFunction::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
33 
34  typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
35  typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
36  typedef typename DiscreteFunctionSpaceType::RangeFieldType RangeFieldType;
37 
38  const int dimDomain = DiscreteFunctionSpaceType::dimDomain;
39  const int dimRange = DiscreteFunctionSpaceType::dimRange;
40 
41  const int size = FourierFunctionSetSize< dimDomain, DiscreteFunctionSpaceType::polynomialOrder-1 >::v;
42 
43  // mass matrix and right hand side for L2 projection
44  FieldMatrix< RangeFieldType, size, size > mat( 0 );
45  FieldMatrix< RangeFieldType, dimRange, size > rhs( 0 );
46 
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 );
51 
52  // Note: The partition is ignored, here. As all degrees of freedom are
53  // global, we need to fill in all DoFs anyway.
54  // As we need to compute the global mass matrix, we only add our
55  // interior contibutions, though.
56  for( const auto entity : elements( v.gridPart(), Partitions::interior ) )
57  {
58  const auto geometry = entity.geometry();
59 
60  uLocal.init( entity );
61 
62  for( const auto qp : Dune::Fem::CachingQuadrature< GridPartType, 0 >( entity, 2*functionSet.order() ) )
63  {
64  // obtain global position and weight
65  const auto x = geometry.global( qp.position() );
66  const RangeFieldType weight = geometry.integrationElement( x ) * qp.weight();
67 
68  // evaluate scalar function set
69  FieldVector< RangeFieldType, size > values;
70  functionSet.evaluateEach( x, [ &values ] ( std::size_t i, const RangeFieldType &v ) { values[ i ] = v; } );
71 
72  // update mass matrix
73  for( int i = 0; i < size; ++i )
74  mat[ i ].axpy( values[ i ]*weight, values );
75 
76  // evaluate u
77  RangeType uValue;
78  uLocal.evaluate( qp, uValue );
79 
80  // update right hand side
81  for( int i = 0; i < dimRange; ++i )
82  rhs[ i ].axpy( uValue[ i ]*weight, values );
83  }
84  }
85 
86  // globally sum up mat and rhs
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() );
97 
98  // solve mat * dofs^T = rhs^T
99  mat.invert();
100  FieldMatrix< RangeFieldType, dimRange, size > dofs( 0 );
101  for( int i = 0; i < dimRange; ++i )
102  mat.umv( rhs[ i ], dofs[ i ] );
103 
104  // copy dofs to dof vector of v
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 ];
108  }
109 
110  } // namespace Fem
111 
112 } // namespace Dune
113 
114 #endif // #ifndef DUNE_FEM_SPACE_FOURIER_INTERPOLATE_HH
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
static constexpr T max(T a)
Definition: utility.hh:77
Definition: space/fourier/functionset.hh:28