dune-fem  2.8-git
discontinuousgalerkin/localinterpolation.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_SPACE_DISCONTINUOUSGALERKIN_LOCALINTERPOLATION_HH
2 #define DUNE_FEM_SPACE_DISCONTINUOUSGALERKIN_LOCALINTERPOLATION_HH
3 
4 // dune-fem includes
5 #include <dune/grid/common/capabilities.hh>
9 
16 namespace Dune
17 {
18 
19  namespace Fem
20  {
21 
22  // DiscontinuousGalerkinLocalInterpolation
23  // ---------------------------------------
24 
28  template< class DiscreteFunctionSpace >
30  {
32 
33  public:
35  typedef typename DiscreteFunctionSpaceType::GridType GridType;
36  typedef typename DiscreteFunctionSpaceType::EntityType EntityType;
37 
38  static const bool isAlwaysAffine = Dune::Capabilities::isCartesian< GridType >::v ||
39  ( Dune::Capabilities::hasSingleGeometryType< GridType >::v && ((Dune::Capabilities::hasSingleGeometryType< GridType >::topologyId >> 1) == 0)) ;
40  // always true for orthonormal spaces
41  //static const bool isAlwaysAffine = true;
42 
43  private:
44  typedef typename DiscreteFunctionSpaceType::RangeType RangeType;
45  typedef typename RangeType::value_type RangeFieldType;
46 
47  typedef typename DiscreteFunctionSpaceType::GridPartType GridPartType;
48 
49  typedef typename DiscreteFunctionSpaceType :: LocalMassMatrixType
50  LocalMassMatrixType;
51 
52  typedef typename LocalMassMatrixType::VolumeQuadratureType QuadratureType;
53 
54  public:
56  : space_( space ),
57  massStorage_( space_.localMassMatrixStorage() )
58  {}
59 
62 
63  void bind( const EntityType& ) {}
64  void unbind() {}
65 
66  ThisType &operator= ( const ThisType &other ) = delete;
67 
68  template< class LocalFunction, class LocalDofVector >
69  void operator () ( const LocalFunction &localFunction, LocalDofVector &dofs ) const
70  {
71  // set all dofs to zero
72  std::fill( dofs.begin(), dofs.end(), typename LocalDofVector::value_type(0) );
73 
74  // get entity and geometry
75  const EntityType &entity = localFunction.entity();
76 
77  if( entity.type().isNone() )
78  {
80  ElementQuadratureType quadrature( entity, localFunction.order() + space_.order( entity ) );
81  bool isAffine = computeInterpolation( entity, quadrature, localFunction, dofs );
82  if( ! isAffine )
83  {
85  AggloMassMatrix massMat( space_, massMatrix().volumeQuadratureOrder( entity ) );
86  // apply inverse of mass matrix
87  auto basisFunctionSet = space_.basisFunctionSet(entity);
88  massMat.applyInverse( entity, basisFunctionSet, dofs );
89  }
90  }
91  else
92  {
93  QuadratureType quadrature( entity, localFunction.order() + space_.order( entity ) );
94  bool isAffine = computeInterpolation( entity, quadrature, localFunction, dofs );
95  if( ! isAffine )
96  {
97  // apply inverse of mass matrix
98  auto basisFunctionSet = space_.basisFunctionSet(entity);
99  massMatrix().applyInverse( entity, basisFunctionSet, dofs );
100  }
101  }
102 
103  }
104 
105  private:
106  template<class QuadImpl, class LocalFunction, class LocalDofVector >
107  bool computeInterpolation( const EntityType& entity,
108  const QuadImpl& quadrature,
110  LocalDofVector &dofs ) const
111  {
112  const int nop = quadrature.nop();
113 
114  auto& values = massStorage_.second;
115  // adjust size of values
116  values.resize( nop );
117 
118  // evaluate local function for all quadrature points
119  localFunction.evaluateQuadrature( quadrature, values );
120 
121  bool isAffine = isAlwaysAffine ;
122  if( ! isAlwaysAffine )
123  {
124  const auto geometry = entity.geometry();
125  isAffine = geometry.affine();
126 
127  if( ! isAffine )
128  {
129  // apply weight
130  for(auto qp : quadrature )
131  values[ qp.index() ] *= qp.weight() * geometry.integrationElement( qp.position() );
132  }
133  }
134 
135  if( isAffine )
136  {
137  // apply weight only
138  for(auto qp : quadrature )
139  values[ qp.index() ] *= qp.weight();
140  }
141 
142  // add values to local function
143  space_.basisFunctionSet(entity).axpy( quadrature, values, dofs );
144 
145  return isAffine;
146  }
147 
148  const LocalMassMatrixType &massMatrix () const
149  {
150  return massStorage_.first;
151  }
152 
153  const DiscreteFunctionSpaceType& space_;
154  typename DiscreteFunctionSpaceType::LocalMassMatrixStorageType& massStorage_;
155  };
156 
157  } // namespace Fem
158 
159 } // namespace Dune
160 
161 #endif // #ifndef DUNE_FEM_SPACE_DISCONTINUOUSGALERKIN_LOCALINTERPOLATION_HH
Definition: bindguard.hh:11
static GridFunctionView< GF > localFunction(const GF &gf)
Definition: gridfunctionview.hh:118
IteratorRange< typename DF::DofIteratorType > dofs(DF &df)
Iterates over all DOFs.
Definition: rangegenerators.hh:76
interface for local functions
Definition: localfunction.hh:77
Local Mass Matrix for arbitrary spaces.
Definition: localmassmatrix.hh:909
quadrature on the codim-0 reference element
Definition: elementquadrature.hh:18
Definition: discontinuousgalerkin/localinterpolation.hh:30
void operator()(const LocalFunction &localFunction, LocalDofVector &dofs) const
Definition: discontinuousgalerkin/localinterpolation.hh:69
void unbind()
Definition: discontinuousgalerkin/localinterpolation.hh:64
void bind(const EntityType &)
Definition: discontinuousgalerkin/localinterpolation.hh:63
DiscontinuousGalerkinLocalInterpolation(const ThisType &other)=default
DiscreteFunctionSpaceType::EntityType EntityType
Definition: discontinuousgalerkin/localinterpolation.hh:36
DiscreteFunctionSpace DiscreteFunctionSpaceType
Definition: discontinuousgalerkin/localinterpolation.hh:34
static const bool isAlwaysAffine
Definition: discontinuousgalerkin/localinterpolation.hh:38
ThisType & operator=(const ThisType &other)=delete
DiscreteFunctionSpaceType::GridType GridType
Definition: discontinuousgalerkin/localinterpolation.hh:35
DiscontinuousGalerkinLocalInterpolation(const DiscreteFunctionSpaceType &space)
Definition: discontinuousgalerkin/localinterpolation.hh:55
DiscontinuousGalerkinLocalInterpolation(ThisType &&other)=default
discrete function space