dune-fem  2.8-git
raviartthomas/localinterpolation.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_SPACE_RAVIARTTHOMAS_LOCALINTERPOLATION_HH
2 #define DUNE_FEM_SPACE_RAVIARTTHOMAS_LOCALINTERPOLATION_HH
3 
4 // alternative interpolation used for testing
5 
6 // C++ includes
7 #include <cassert>
8 #include <vector>
9 #include <utility>
10 
11 // dune-common includes
12 #include <dune/common/fvector.hh>
13 #include <dune/common/typetraits.hh>
14 
15 // dune-geometry includes
16 #include <dune/geometry/type.hh>
17 #include <dune/geometry/referenceelements.hh>
18 #include <dune/geometry/quadraturerules.hh>
19 
20 // dune-fem includes
22 
23 
24 namespace Dune
25 {
26 
27  namespace Fem
28  {
29 
30  namespace Impl
31  {
32 
33  // forward declarations
34  // --------------------
35 
36  template< unsigned int, class, class, int, int > struct RaviartThomasLocalFiniteElement;
37 
38 
39  // RaviartThomasLocalInterpolationBasis
40  // ------------------------------------
41 
42  /*
43  * These are mostly copies from the interpolation implementations in dune-localfunctions
44  * (dune/localfunctions/raviartthomas/raviartthomas * /raviartthomas * localinterpolation.hh)
45  */
46 
47  template< class LocalFiniteElement >
48  struct RaviartThomasLocalInterpolationBasis
49  {
50  static_assert( AlwaysFalse< LocalFiniteElement >::value,
51  "`RaviartThomasLocalInterpolationBasis` not implemented for your choice." );
52  };
53 
54 
55  // 0th order
56  template< unsigned int id, class Domain, class Range, int dim >
57  struct RaviartThomasLocalInterpolationBasis< RaviartThomasLocalFiniteElement< id, Domain, Range, dim, 0 > >
58  {
59  using DomainType = FieldVector< Domain, dim >;
60  using FaceDomainType = FieldVector< Domain, dim-1 >;
61  using RangeType = FieldVector< Range, dim >;
62  using RangeFieldType = Range;
63 
64  RaviartThomasLocalInterpolationBasis () = default;
65  explicit RaviartThomasLocalInterpolationBasis ( unsigned int orientations ) : orientations_( orientations ) {}
66 
67  void trace ( int facet, const FaceDomainType& x, std::vector< std::pair< int, RangeFieldType > >& basis ) const
68  {
69  assert( basis.size() >= size( 1 ) );
70  basis[ 0 ] = std::make_pair( facet, sign( facet ) );
71  }
72 
73  void interior ( const DomainType& x, std::vector< std::pair< int, RangeType > >& basis ) const {}
74 
75  constexpr auto size ( int codim ) const -> std::size_t { assert( codim < 2 ); return (codim == 0) ? 0 : 1; }
76  constexpr int order ( int codim ) const { assert( codim < 2 ); return (codim == 0) ? 0 : 1; }
77 
78  private:
79  auto sign ( int facet ) const -> RangeFieldType { return (orientations_ & (1u << facet)) ? -1.0 : 1.0; }
80 
81  unsigned int orientations_ = 0;
82  };
83 
84  // 1st order, 2d Simplex
85  template< class Domain, class Range >
86  struct RaviartThomasLocalInterpolationBasis< RaviartThomasLocalFiniteElement< Dune::Impl::SimplexTopology< 2 >::type::id, Domain, Range, 2, 1 > >
87  {
88  using DomainType = FieldVector< Domain, 2 >;
89  using FaceDomainType = FieldVector< Domain, 1 >;
90  using RangeType = FieldVector< Range, 2 >;
91  using RangeFieldType = Range;
92 
93  RaviartThomasLocalInterpolationBasis () = default;
94  explicit RaviartThomasLocalInterpolationBasis ( unsigned int orientations ) : orientations_( orientations ) {}
95 
96  void trace ( int facet, const FaceDomainType& x, std::vector< std::pair< int, RangeFieldType > >& basis ) const
97  {
98  assert( basis.size() >= size( 1 ) );
99  const RangeFieldType temp = (facet==1) ? 1.0 - 2.0*x : 2.0*x - 1.0;
100  basis[ 0 ] = std::make_pair( facet, sign( facet ) );
101  basis[ 1 ] = std::make_pair( facet+3, temp );
102  }
103 
104  void interior ( const DomainType& x, std::vector< std::pair< int, RangeType > >& basis ) const
105  {
106  assert( basis.size() >= size( 0 ) );
107  basis[ 0 ] = std::make_pair( 6, RangeType{ 1.0, 0.0 } );
108  basis[ 1 ] = std::make_pair( 7, RangeType{ 0.0, 1.0 } );
109  }
110 
111  constexpr auto size ( int codim ) const -> std::size_t { assert( codim < 2 ); return (codim == 0) ? 2 : 2; }
112  constexpr int order ( int codim ) const { assert( codim < 2 ); return (codim == 0) ? 8 : 4; }
113 
114 
115  private:
116  auto sign ( int facet ) const -> RangeFieldType { return (orientations_ & (1u << facet)) ? -1.0 : 1.0; }
117 
118  unsigned int orientations_ = 0;
119  };
120 
121  // 1st order, 2d Cube
122  template< class Domain, class Range >
123  struct RaviartThomasLocalInterpolationBasis< RaviartThomasLocalFiniteElement< Dune::GeometryTypes::cube( 2 ).id(), Domain, Range, 2, 1 > >
124  {
125  using DomainType = FieldVector< Domain, 2 >;
126  using FaceDomainType = FieldVector< Domain, 1 >;
127  using RangeType = FieldVector< Range, 2 >;
128  using RangeFieldType = Range;
129 
130  RaviartThomasLocalInterpolationBasis () = default;
131  explicit RaviartThomasLocalInterpolationBasis ( unsigned int orientations ) : orientations_( orientations ) {}
132 
133  void trace ( int facet, const FaceDomainType& x, std::vector< std::pair< int, RangeFieldType > >& basis ) const
134  {
135  assert( basis.size() >= size( 1 ) );
136  const RangeFieldType temp = (facet > 0 && facet < 3 ) ? 1.0 - 2.0*x : 2.0*x - 1.0;
137  basis[ 0 ] = std::make_pair( 2*facet , sign( facet ) );
138  basis[ 1 ] = std::make_pair( 2*facet+1, temp );
139  }
140 
141  void interior ( const DomainType& x, std::vector< std::pair< int, RangeType > >& basis ) const
142  {
143  assert( basis.size() >= size( 0 ) );
144  basis[ 0 ] = std::make_pair( 8, RangeType{ 1.0, 0.0 } );
145  basis[ 1 ] = std::make_pair( 9, RangeType{ 0.0, 1.0 } );
146  basis[ 2 ] = std::make_pair( 10, RangeType{ x[1], 0.0 } );
147  basis[ 3 ] = std::make_pair( 11, RangeType{ 0.0, x[0] } );
148  }
149 
150  constexpr auto size ( int codim ) const -> std::size_t { assert( codim < 2 ); return (codim == 0) ? 4 : 2; }
151  constexpr int order ( int codim ) const { assert( codim < 2 ); return (codim == 0) ? 3 : 3; }
152 
153  private:
154  auto sign ( int facet ) const -> RangeFieldType { return (orientations_ & (1u << facet)) ? -1.0 : 1.0; }
155 
156  unsigned int orientations_ = 0;
157  };
158 
159  // 1st order, 3d Cube
160  template< class Domain, class Range >
161  struct RaviartThomasLocalInterpolationBasis< RaviartThomasLocalFiniteElement< Dune::GeometryTypes::cube( 3 ).id(), Domain, Range, 3, 1 > >
162  {
163  using DomainType = FieldVector< Domain, 3 >;
164  using FaceDomainType = FieldVector< Domain, 2 >;
165  using RangeType = FieldVector< Range, 3 >;
166  using RangeFieldType = Range;
167 
168  RaviartThomasLocalInterpolationBasis () = default;
169  explicit RaviartThomasLocalInterpolationBasis ( unsigned int orientations ) : orientations_( orientations ) {}
170 
171  void trace ( int facet, const FaceDomainType& x, std::vector< std::pair< int, RangeFieldType > >& basis ) const
172  {
173  assert( basis.size() >= size( 1 ) );
174 
175  const RangeFieldType tempX = 2.0*x[0] - 1.0;
176  const RangeFieldType tempY = 2.0*x[1] - 1.0;
177 
178  basis[ 0 ] = std::make_pair( facet, sign( facet ) );
179 
180  switch( facet )
181  {
182  case 0:
183  case 5:
184  basis[ 1 ] = std::make_pair( facet+ 6, tempX );
185  basis[ 2 ] = std::make_pair( facet+12, tempY );
186  basis[ 3 ] = std::make_pair( facet+18, tempX*tempY );
187  break;
188  case 1:
189  case 4:
190  basis[ 1 ] = std::make_pair( facet+ 6, -tempX );
191  basis[ 2 ] = std::make_pair( facet+12, -tempY );
192  basis[ 3 ] = std::make_pair( facet+18, -tempX*tempY );
193  break;
194  case 2:
195  basis[ 1 ] = std::make_pair( facet+ 6, -tempX );
196  basis[ 2 ] = std::make_pair( facet+12, tempY );
197  basis[ 3 ] = std::make_pair( facet+18, -tempX*tempY );
198  break;
199  case 3:
200  basis[ 1 ] = std::make_pair( facet+ 6, tempX );
201  basis[ 2 ] = std::make_pair( facet+12, -tempY );
202  basis[ 3 ] = std::make_pair( facet+18, tempX*tempY );
203  break;
204  }
205  }
206 
207  void interior ( const DomainType& x, std::vector< std::pair< int, RangeType > >& basis ) const
208  {
209  assert( basis.size() >= size( 0 ) );
210  basis[ 0 ] = std::make_pair( 24, RangeType{ 1.0, 0.0, 0.0 } );
211  basis[ 1 ] = std::make_pair( 25, RangeType{ 0.0, 1.0, 0.0 } );
212  basis[ 2 ] = std::make_pair( 26, RangeType{ 0.0, 0.0, 1.0 } );
213  basis[ 3 ] = std::make_pair( 27, RangeType{ x[1], 0.0, 0.0 } );
214  basis[ 4 ] = std::make_pair( 28, RangeType{ x[2], 0.0, 0.0 } );
215  basis[ 5 ] = std::make_pair( 29, RangeType{ 0.0, x[0], 0.0 } );
216  basis[ 6 ] = std::make_pair( 30, RangeType{ 0.0, x[2], 0.0 } );
217  basis[ 7 ] = std::make_pair( 31, RangeType{ 0.0, 0.0, x[0] } );
218  basis[ 8 ] = std::make_pair( 32, RangeType{ 0.0, 0.0, x[1] } );
219  basis[ 9 ] = std::make_pair( 33, RangeType{ x[1]*x[2], 0.0, 0.0 } );
220  basis[ 10 ] = std::make_pair( 34, RangeType{ 0.0, x[0]*x[2], 0.0 } );
221  basis[ 11 ] = std::make_pair( 35, RangeType{ 0.0, 0.0, x[0]*x[1] } );
222  }
223 
224  constexpr auto size ( int codim ) const -> std::size_t { assert( codim < 2 ); return (codim == 0) ? 12 : 4; }
225  constexpr int order ( int codim ) const { assert( codim < 2 ); return (codim == 0) ? 3 : 3; }
226 
227  private:
228  auto sign ( int facet ) const -> RangeFieldType { return (orientations_ & (1u << facet)) ? -1.0 : 1.0; }
229 
230  unsigned int orientations_ = 0;
231  };
232 
233  }
234 
235  template< class GridPart, class LocalFiniteElement, int dimRange = GridPart::dimension >
237  {
238  public:
239  using GridPartType = GridPart;
240  using LocalFiniteElementType = LocalFiniteElement;
241 
242  using EntityType = typename GridPartType::template Codim< 0 >::EntityType;
243 
244  protected:
245  using LocalInterpolationBasisType = Impl::RaviartThomasLocalInterpolationBasis< LocalFiniteElementType >;
246 
247  using Geometry = typename EntityType::Geometry;
248  using ReferenceElementType = ReferenceElement< Geometry >;
249 
252 
253  using RangeType = typename LocalInterpolationBasisType::RangeType;
254  using RangeFieldType = typename LocalInterpolationBasisType::RangeFieldType;
255 
256  using VolumeQuadratures = QuadratureRules< typename Geometry::ctype, ReferenceElementType::dimension >;
257  using FaceQuadratures = QuadratureRules< RangeFieldType, ReferenceElementType::dimension-1 >;
258 
259  public:
260  explicit RaviartThomasLocalInterpolation ( const EntityType& entity )
261  : RaviartThomasLocalInterpolation( entity, 0, -1 )
262  {}
263 
264  RaviartThomasLocalInterpolation ( const EntityType& entity, unsigned int orientations, int order = -1 )
265  : geometry_( entity.geometry() ),
266  refElement_( Dune::referenceElement( geometry_ ) ),
267  localBasis_( orientations ),
268  order_( order )
269  {}
270 
271  template< class LocalFunction, class LocalDofVector >
272  void interior ( const LocalFunction& lf, LocalDofVector& dofs ) const
273  {
274  if ( !hasInterior() )
275  return;
276 
277  assert( dofs.size() >= localBasis().size( 0 ) + referenceElement().size( 1 ) * localBasis().size( 1 ) );
278  interior_.resize( localBasis().size( 0 ) );
279 
280  for( const auto& qp : getQuadrature( (order_ == -1) ? localBasis().order( 0 ) : order_ ) )
281  {
282  auto invPiola = inverseTransformation( qp.position() );
283 
284  localBasis().interior( qp.position(), interior_ );
285  auto value = invPiola.apply( lf( qp.position() ) );
286 
287  for ( const auto& base : interior_ )
288  dofs[ base.first ] += qp.weight() * (value * base.second);
289  }
290  }
291 
292  template< class LocalFunction, class LocalDofVector >
293  void trace ( int facet, const LocalFunction& lf, LocalDofVector& dofs ) const
294  {
295  assert( dofs.size() >= localBasis().size( 0 ) + referenceElement().size( 1 ) * localBasis().size( 1 ) );
296  traces_.resize( localBasis().size( 1 ) );
297 
298  auto embedding = referenceElement().template geometry< 1 >( facet );
299  auto normal = referenceElement().integrationOuterNormal( facet );
300 
301  for ( const auto& qp : getQuadrature( facet, (order_ == -1) ? localBasis().order( 1 ) : order_ ) )
302  {
303  auto p = embedding.global( qp.position() );
304  auto invPiola = inverseTransformation( p );
305 
306  localBasis().trace( facet, qp.position(), traces_ );
307  auto value = invPiola.apply( lf( p ) );
308 
309  for ( const auto& base : traces_ )
310  dofs[ base.first ] += qp.weight() * (value * normal) * base.second;
311  }
312  }
313 
314  template< class LocalFunction, class LocalDofVector >
315  void interiorTrace ( int facet, const LocalFunction& lf, LocalDofVector& dofs ) const
316  {
317  if ( !hasInterior() )
318  return;
319 
320  assert( dofs.size() >= localBasis().size( 0 ) + referenceElement().size( 1 ) * localBasis().size( 1 ) );
321  interior_.resize( localBasis().size( 0 ) );
322 
323  auto embedding = referenceElement().template geometry< 1 >( facet );
324 
325  for( const auto& qp : getQuadrature( facet, (order_ == -1) ? localBasis().order( 0 ) : order_ ) )
326  {
327  auto p = embedding.global( qp.position() );
328  auto invPiola = inverseTransformation( p );
329 
330  localBasis().interior( p, interior_ );
331  auto value = invPiola.apply( lf( p ) );
332 
333  for ( const auto& base : interior_ )
334  dofs[ base.first ] += qp.weight() * (value * base.second);
335  }
336  }
337 
338  template< class LocalFunction, class LocalDofVector >
339  void operator() ( const LocalFunction& lf, LocalDofVector& dofs ) const
340  {
341  for ( int facet : range( referenceElement().size( 1 ) ) )
342  trace( facet, lf, dofs );
343  interior( lf, dofs );
344  }
345 
346  bool hasInterior () const { return localBasis().size( 0 ) > 0; }
347 
348  protected:
349  auto geometry () const -> const Geometry& { return geometry_; }
350  auto referenceElement () const -> const ReferenceElementType& { return refElement_; }
351  auto localBasis () const -> const LocalInterpolationBasisType& { return localBasis_; }
352 
353  template< class Point >
354  auto transformation ( const Point& p ) const { return TransformationType( geometry(), p ); }
355 
356  template< class Point >
357  auto inverseTransformation ( const Point& p ) const { return InverseTransformationType( geometry(), p ); }
358 
359  auto getQuadrature ( int order ) const { return VolumeQuadratures::rule( referenceElement().type(), order ); }
360  auto getQuadrature ( int facet, int order ) const { return FaceQuadratures::rule( referenceElement().type( facet, 1 ), order ); }
361 
362  private:
363  Geometry geometry_;
364  ReferenceElementType refElement_;
365  LocalInterpolationBasisType localBasis_;
366  const int order_;
367 
368  mutable std::vector< std::pair< int, RangeFieldType > > traces_;
369  mutable std::vector< std::pair< int, RangeType > > interior_;
370  };
371 
372  } // namespace Fem
373 
374 } // namespace Dune
375 
376 #endif // #ifndef DUNE_FEM_SPACE_RAVIARTTHOMAS_LOCALINTERPOLATION_HH
Definition: bindguard.hh:11
IteratorRange< typename DF::DofIteratorType > dofs(DF &df)
Iterates over all DOFs.
Definition: rangegenerators.hh:76
interface for local functions
Definition: localfunction.hh:77
Definition: piolatransformation.hh:48
InversePiolaTransformation< Geometry, dimRange > InverseTransformationType
Definition: piolatransformation.hh:60
Definition: raviartthomas/localinterpolation.hh:237
auto inverseTransformation(const Point &p) const
Definition: raviartthomas/localinterpolation.hh:357
typename TransformationType::InverseTransformationType InverseTransformationType
Definition: raviartthomas/localinterpolation.hh:251
auto getQuadrature(int order) const
Definition: raviartthomas/localinterpolation.hh:359
typename EntityType::Geometry Geometry
Definition: raviartthomas/localinterpolation.hh:247
QuadratureRules< RangeFieldType, ReferenceElementType::dimension-1 > FaceQuadratures
Definition: raviartthomas/localinterpolation.hh:257
auto getQuadrature(int facet, int order) const
Definition: raviartthomas/localinterpolation.hh:360
typename GridPartType::template Codim< 0 >::EntityType EntityType
Definition: raviartthomas/localinterpolation.hh:242
ReferenceElement< Geometry > ReferenceElementType
Definition: raviartthomas/localinterpolation.hh:248
void trace(int facet, const LocalFunction &lf, LocalDofVector &dofs) const
Definition: raviartthomas/localinterpolation.hh:293
QuadratureRules< typename Geometry::ctype, ReferenceElementType::dimension > VolumeQuadratures
Definition: raviartthomas/localinterpolation.hh:256
void interiorTrace(int facet, const LocalFunction &lf, LocalDofVector &dofs) const
Definition: raviartthomas/localinterpolation.hh:315
auto geometry() const -> const Geometry &
Definition: raviartthomas/localinterpolation.hh:349
Impl::RaviartThomasLocalInterpolationBasis< LocalFiniteElementType > LocalInterpolationBasisType
Definition: raviartthomas/localinterpolation.hh:245
auto transformation(const Point &p) const
Definition: raviartthomas/localinterpolation.hh:354
PiolaTransformation< Geometry, dimRange > TransformationType
Definition: raviartthomas/localinterpolation.hh:250
GridPart GridPartType
Definition: raviartthomas/localinterpolation.hh:239
typename LocalInterpolationBasisType::RangeFieldType RangeFieldType
Definition: raviartthomas/localinterpolation.hh:254
auto localBasis() const -> const LocalInterpolationBasisType &
Definition: raviartthomas/localinterpolation.hh:351
typename LocalInterpolationBasisType::RangeType RangeType
Definition: raviartthomas/localinterpolation.hh:253
LocalFiniteElement LocalFiniteElementType
Definition: raviartthomas/localinterpolation.hh:240
void interior(const LocalFunction &lf, LocalDofVector &dofs) const
Definition: raviartthomas/localinterpolation.hh:272
void operator()(const LocalFunction &lf, LocalDofVector &dofs) const
Definition: raviartthomas/localinterpolation.hh:339
RaviartThomasLocalInterpolation(const EntityType &entity)
Definition: raviartthomas/localinterpolation.hh:260
auto referenceElement() const -> const ReferenceElementType &
Definition: raviartthomas/localinterpolation.hh:350
RaviartThomasLocalInterpolation(const EntityType &entity, unsigned int orientations, int order=-1)
Definition: raviartthomas/localinterpolation.hh:264
bool hasInterior() const
Definition: raviartthomas/localinterpolation.hh:346