DUNE-FEM (unstable)

localinterpolation.hh
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
21#include <dune/fem/space/basisfunctionset/piolatransformation.hh>
22
23
24namespace 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::GeometryTypes::simplex(2).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 >
236 class RaviartThomasLocalInterpolation
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
250 using TransformationType = PiolaTransformation< Geometry, dimRange >;
251 using InverseTransformationType = typename TransformationType::InverseTransformationType;
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
Dune::QuadratureRules::rule
static const QuadratureRule & rule(const GeometryType &t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:326
fvector.hh
Implements a vector constructed from a given type representing a field and a compile-time given size.
Dune::referenceElement
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
Dune::GeometryTypes::cube
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:462
Dune::GeometryTypes::simplex
constexpr GeometryType simplex(unsigned int dim)
Returns a GeometryType representing a simplex of dimension dim.
Definition: type.hh:453
Dune::Partitions::interior
constexpr Interior interior
PartitionSet for the interior partition.
Definition: partitionset.hh:271
Dune
Dune namespace.
Definition: alignedallocator.hh:13
Dune::size
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
Dune::sign
int sign(const T &val)
Return the sign of the value.
Definition: math.hh:180
type.hh
A unique label for each type of element that can occur in a grid.
typetraits.hh
Traits for type conversions and type information.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 27, 22:29, 2024)