1#ifndef DUNE_FEM_SPACE_RAVIARTTHOMAS_LOCALINTERPOLATION_HH
2#define DUNE_FEM_SPACE_RAVIARTTHOMAS_LOCALINTERPOLATION_HH
12#include <dune/common/fvector.hh>
13#include <dune/common/typetraits.hh>
16#include <dune/geometry/type.hh>
17#include <dune/geometry/referenceelements.hh>
18#include <dune/geometry/quadraturerules.hh>
36 template<
unsigned int,
class,
class,
int,
int >
struct RaviartThomasLocalFiniteElement;
47 template<
class LocalFiniteElement >
48 struct RaviartThomasLocalInterpolationBasis
50 static_assert( AlwaysFalse< LocalFiniteElement >::value,
51 "`RaviartThomasLocalInterpolationBasis` not implemented for your choice." );
56 template<
unsigned int id,
class Domain,
class Range,
int dim >
57 struct RaviartThomasLocalInterpolationBasis< RaviartThomasLocalFiniteElement< id, Domain, Range, dim, 0 > >
59 using DomainType = FieldVector< Domain, dim >;
60 using FaceDomainType = FieldVector< Domain, dim-1 >;
61 using RangeType = FieldVector< Range, dim >;
62 using RangeFieldType = Range;
64 RaviartThomasLocalInterpolationBasis () =
default;
65 explicit RaviartThomasLocalInterpolationBasis (
unsigned int orientations ) : orientations_( orientations ) {}
67 void trace (
int facet,
const FaceDomainType& x, std::vector< std::pair< int, RangeFieldType > >& basis )
const
69 assert( basis.size() >= size( 1 ) );
70 basis[ 0 ] = std::make_pair( facet, sign( facet ) );
73 void interior (
const DomainType& x, std::vector< std::pair< int, RangeType > >& basis )
const {}
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; }
79 auto sign (
int facet )
const -> RangeFieldType {
return (orientations_ & (1u << facet)) ? -1.0 : 1.0; }
81 unsigned int orientations_ = 0;
85 template<
class Domain,
class Range >
86 struct RaviartThomasLocalInterpolationBasis< RaviartThomasLocalFiniteElement<
Dune::Impl::SimplexTopology< 2 >::type::id, Domain, Range, 2, 1 > >
88 using DomainType = FieldVector< Domain, 2 >;
89 using FaceDomainType = FieldVector< Domain, 1 >;
90 using RangeType = FieldVector< Range, 2 >;
91 using RangeFieldType = Range;
93 RaviartThomasLocalInterpolationBasis () =
default;
94 explicit RaviartThomasLocalInterpolationBasis (
unsigned int orientations ) : orientations_( orientations ) {}
96 void trace (
int facet,
const FaceDomainType& x, std::vector< std::pair< int, RangeFieldType > >& basis )
const
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 );
104 void interior (
const DomainType& x, std::vector< std::pair< int, RangeType > >& basis )
const
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 } );
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; }
116 auto sign (
int facet )
const -> RangeFieldType {
return (orientations_ & (1u << facet)) ? -1.0 : 1.0; }
118 unsigned int orientations_ = 0;
122 template<
class Domain,
class Range >
123 struct RaviartThomasLocalInterpolationBasis< RaviartThomasLocalFiniteElement<
Dune::GeometryTypes::cube( 2 ).id(), Domain, Range, 2, 1 > >
125 using DomainType = FieldVector< Domain, 2 >;
126 using FaceDomainType = FieldVector< Domain, 1 >;
127 using RangeType = FieldVector< Range, 2 >;
128 using RangeFieldType = Range;
130 RaviartThomasLocalInterpolationBasis () =
default;
131 explicit RaviartThomasLocalInterpolationBasis (
unsigned int orientations ) : orientations_( orientations ) {}
133 void trace (
int facet,
const FaceDomainType& x, std::vector< std::pair< int, RangeFieldType > >& basis )
const
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 );
141 void interior (
const DomainType& x, std::vector< std::pair< int, RangeType > >& basis )
const
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] } );
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; }
154 auto sign (
int facet )
const -> RangeFieldType {
return (orientations_ & (1u << facet)) ? -1.0 : 1.0; }
156 unsigned int orientations_ = 0;
160 template<
class Domain,
class Range >
161 struct RaviartThomasLocalInterpolationBasis< RaviartThomasLocalFiniteElement<
Dune::GeometryTypes::cube( 3 ).id(), Domain, Range, 3, 1 > >
163 using DomainType = FieldVector< Domain, 3 >;
164 using FaceDomainType = FieldVector< Domain, 2 >;
165 using RangeType = FieldVector< Range, 3 >;
166 using RangeFieldType = Range;
168 RaviartThomasLocalInterpolationBasis () =
default;
169 explicit RaviartThomasLocalInterpolationBasis (
unsigned int orientations ) : orientations_( orientations ) {}
171 void trace (
int facet,
const FaceDomainType& x, std::vector< std::pair< int, RangeFieldType > >& basis )
const
173 assert( basis.size() >= size( 1 ) );
175 const RangeFieldType tempX = 2.0*x[0] - 1.0;
176 const RangeFieldType tempY = 2.0*x[1] - 1.0;
178 basis[ 0 ] = std::make_pair( facet, sign( facet ) );
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 );
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 );
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 );
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 );
207 void interior (
const DomainType& x, std::vector< std::pair< int, RangeType > >& basis )
const
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] } );
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; }
228 auto sign (
int facet )
const -> RangeFieldType {
return (orientations_ & (1u << facet)) ? -1.0 : 1.0; }
230 unsigned int orientations_ = 0;
235 template<
class Gr
idPart,
class LocalFiniteElement,
int dimRange = Gr
idPart::dimension >
242 using EntityType =
typename GridPartType::template Codim< 0 >::EntityType;
253 using RangeType =
typename LocalInterpolationBasisType::RangeType;
256 using VolumeQuadratures = QuadratureRules< typename Geometry::ctype, ReferenceElementType::dimension >;
267 localBasis_( orientations ),
271 template<
class LocalFunction,
class LocalDofVector >
284 localBasis().interior( qp.position(), interior_ );
285 auto value = invPiola.apply( lf( qp.position() ) );
287 for (
const auto& base : interior_ )
288 dofs[ base.first ] += qp.weight() * (value * base.second);
292 template<
class LocalFunction,
class LocalDofVector >
303 auto p = embedding.global( qp.position() );
306 localBasis().trace( facet, qp.position(), traces_ );
307 auto value = invPiola.apply( lf( p ) );
309 for (
const auto& base : traces_ )
310 dofs[ base.first ] += qp.weight() * (value * normal) * base.second;
314 template<
class LocalFunction,
class LocalDofVector >
327 auto p = embedding.global( qp.position() );
331 auto value = invPiola.apply( lf( p ) );
333 for (
const auto& base : interior_ )
334 dofs[ base.first ] += qp.weight() * (value * base.second);
338 template<
class LocalFunction,
class LocalDofVector >
353 template<
class Po
int >
356 template<
class Po
int >
368 mutable std::vector< std::pair< int, RangeFieldType > > traces_;
369 mutable std::vector< std::pair< int, RangeType > > interior_;
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