dune-fem 2.8.0
Loading...
Searching...
No Matches
p1bubble.hh
Go to the documentation of this file.
1#ifndef SPACE_P1BUBBLE_HH
2#define SPACE_P1BUBBLE_HH
3
4#include <vector>
5#include <dune/common/exceptions.hh>
6#include <dune/geometry/referenceelements.hh>
7#include <dune/grid/common/gridenums.hh>
8
20
22
23namespace Dune
24{
25
26 namespace Fem
27 {
28
29 template< class FunctionSpace >
31 {
36
38 : points_( dimDomain + 2, DomainType( 0.0 ) )
39 {
40 for( int i = 0; i < dimDomain; ++i )
41 points_[ i + 1 ][ i ] = 1.0;
42
43 points_[ dimDomain +1 ] = DomainType( 1.0 / ( dimDomain + 1.0 ) );
44 }
45
48
50 template< class LocalFunction, class LocalDofVector >
51 void operator() ( const LocalFunction &lf, LocalDofVector &ldv ) const
53 {
54 int k = 0;
55 for( const DomainType &x : points_ )
56 {
57 RangeType phi;
58 lf.evaluate( x, phi );
59 for( int i = 0; i < dimRange; ++i )
60 ldv[ k++ ] = phi[ i ];
61 }
62 }
63
64 template <class Entity>
65 void bind( const Entity & ) {}
66 void unbind() {}
67
68 private:
69 std::vector< DomainType > points_;
70 };
71
72
73 template< int dim >
75 {
77 BubbleElementLocalKeyMap ( int vertices )
78 {
79 for( int i = 0; i < vertices; ++i )
80 map_.emplace_back( i, dim, 0 );
81 map_.emplace_back( 0, 0, 0 );
82 }
84
85 std::size_t size() const { return map_.size(); }
86
87 LocalKey& localKey ( std::size_t i ) { return map_[ i ]; }
88 const LocalKey& localKey ( std::size_t i ) const { return map_[ i ]; }
89
90 private:
91 std::vector< LocalKey > map_;
92 };
93
95 {
96 // return the shape functions for a given reference element. If this
97 // is not possible an empty DofMapperCode is returned.
98 template< class RefElement,
99 std::enable_if_t< std::is_same< std::decay_t< decltype( std::declval< const RefElement & >().size( 0 ) ) >, int >::value, int > = 0,
100 std::enable_if_t< std::is_same< std::decay_t< decltype( std::declval< const RefElement & >().type( 0, 0 ) ) >, GeometryType >::value, int > = 0 >
101 DofMapperCode operator() ( const RefElement &refElement ) const
102 {
103 static const int dim = RefElement::dimension;
104 if( refElement.type().isSimplex() )
105 return compile( refElement, BubbleElementLocalKeyMap< dim >(dim+1) );
106 if( refElement.type().isCube() && refElement.type().dim() == 2)
107 return compile( refElement, BubbleElementLocalKeyMap< dim >(pow(2,dim)) );
108 else
109 return DofMapperCode();
110 }
111 };
112
113 // SimplexBubbleElementShapeFunctionSet
114 // ----------------------------------
115
116 template< class FunctionSpace >
118 {
120 public:
121 static const int dimDomain = FunctionSpace :: dimDomain;
122 static const int polynomialOrder = dimDomain + 1;
123 static const int numShapeFunctions = dimDomain + 2;
124
126 typedef typename FunctionSpace :: DomainType DomainType;
127 typedef typename FunctionSpace :: RangeType RangeType;
128 static_assert( RangeType::dimension == 1, "This is a scalar shapefunction set" );
129 typedef typename FunctionSpace :: JacobianRangeType JacobianRangeType;
130 typedef typename FunctionSpace :: HessianRangeType HessianRangeType;
131
133
134 template<class GeometryType >
135 SimplexBubbleElementShapeFunctionSet ( const GeometryType& gt )
136 {
137 if( !gt.isSimplex() )
138 DUNE_THROW( NotImplemented, "Wrong geometry type for Cube2DBubblelementShapeFunctionSet." );
139 }
140
141 template< class Point, class Functor >
142 void evaluateEach ( const Point &x, Functor functor ) const
143 {
144 DomainType xRef = coordinate( x );
145 RangeType phi(1), phi0(1);
146 for( int i=0; i< dimDomain; ++i )
147 {
148 functor( i+1, RangeType( xRef[ i ] ) );
149 phi0[ 0 ] -= xRef[ i ];
150 phi[ 0 ] *= xRef[ i ] ;
151 }
152
153 phi[ 0 ] *= phi0[ 0 ] / std::pow( ( dimDomain + 1.0 ), dimDomain + 1.0 );
154 functor( 0, phi0 );
155 functor( dimDomain +1, phi );
156 }
157
158 template< class Point, class Functor >
159 void jacobianEach ( const Point &x, Functor functor ) const
160 {
161 DomainType xRef = coordinate( x );
162
163 JacobianRangeType jac(0), jac0( -1 );
164 RangeType phi0( 1 );
165
166 functor( 0, jac0 );
167
168 for( int i=0; i< dimDomain; ++i )
169 {
170 phi0[ 0 ] -= xRef[ i ];
171
172 for( int j=1; j < dimDomain; ++j )
173 jac0[ 0 ][ (i+j)%dimDomain ] *= xRef[ i ];
174
175 jac[ 0 ][ i ] = 1;
176 functor( i+1, jac );
177 jac[ 0 ][ i ] = 0;
178 }
179
180 for( int i=0; i< dimDomain; ++i )
181 jac0[ 0 ][ i ] *= -(phi0[ 0 ] - xRef[ i ]);
182 jac0[ 0 ] *= 1.0 / std::pow( dimDomain + 1.0, dimDomain + 1.0 );
183 functor( dimDomain +1, jac0 );
184 }
185
186 template< class Point, class Functor >
187 void hessianEach ( const Point &x, Functor functor ) const
188 {
189 DUNE_THROW( NotImplemented, "NotImplemented" );
190 DomainType xRef = coordinate( x );
192 functor( 0, hes );
193 functor( 1, hes );
194 functor( 2, hes );
195 functor( 3, hes );
196 }
197
198 int order () const { return dimDomain + 1; }
199
200 std::size_t size () const { return dimDomain +2; }
201 };
202
203 // 2d cube element taken from
204 // http://arxiv.org/pdf/1507.04417v1.pdf
205 template< class FunctionSpace >
207 {
209 public:
210 static const int dimDomain = FunctionSpace :: dimDomain;
211 static const int polynomialOrder = dimDomain + 1;
212 static const int numShapeFunctions = dimDomain + 2;
213
215 typedef typename FunctionSpace :: DomainType DomainType;
216 typedef typename FunctionSpace :: RangeType RangeType;
217 static_assert( RangeType::dimension == 1, "This is a scalar shapefunction set" );
218 typedef typename FunctionSpace :: JacobianRangeType JacobianRangeType;
219 typedef typename FunctionSpace :: HessianRangeType HessianRangeType;
220
223
224 template<class GeometryType >
225 Cube2DBubbleElementShapeFunctionSet ( const GeometryType& gt )
226 {
227 if( !gt.isCube() )
228 DUNE_THROW( NotImplemented, "Wrong geometry type for Cube2DBubblelementShapeFunctionSet." );
229 if( gt.dim() != 2 )
230 DUNE_THROW( NotImplemented, "2DCubeBubbleElementShapeFunctionSet only implemented for dimension = 2." );
231 }
232
233 template< class Point, class Functor >
234 void evaluateEach ( const Point &x, Functor functor ) const
236 {
237 DomainType xRef = coordinate( x );
238 RangeType phi(1), phi0(1);
239 // (1-x)(1-y) -> grad = ( (y-1),(x-1) )
240 phi0[0] = (1.-xRef[0])*(1.-xRef[1]);
241 functor( 0, phi0 );
242 // x(1-y) -> grad = ( (1-y),-x )
243 phi[0] = xRef[0]*(1.-xRef[1]);
244 functor( 1, phi );
245 // (1-x)y -> grad = ( (-y,(1-x) )
246 phi[0] = (1.-xRef[0])*xRef[1];
247 functor( 2, phi );
248 // xy -> grad = ( (y,x) )
249 phi[0] = xRef[0]*xRef[1];
250 functor( 3, phi );
251 // 64 xy phi_0(x,y)^2 -> grad = 128 phi_0 xy grad_0 +
252 // 64 phi_0^2 (y,x)
253 phi[0] = 64.*phi0[0]*phi0[0]*xRef[0]*xRef[1];
254 functor( 4, phi );
255 }
256
257 template< class Point, class Functor >
258 void jacobianEach ( const Point &x, Functor functor ) const
259 {
260 DomainType xRef = coordinate( x );
261
262 JacobianRangeType jac, jac0;
263 RangeType phi0;
264 phi0[0] = (1.-xRef[0])*(1.-xRef[1]);
265 jac0[0] = {xRef[1]-1,xRef[0]-1};
266
267 functor( 0, jac0 );
268 jac[0] = {1-xRef[1],xRef[0]};
269 functor( 1, jac );
270 jac[0] = {xRef[1],1-xRef[0]};
271 functor( 2, jac );
272 jac[0] = {xRef[1],xRef[0]};
273 functor( 3, jac );
274 // 64 xy phi_0(x,y)^2 -> grad = 128 phi_0 xy grad_0 +
275 // 64 phi_0^2 (y,x)
276 jac0[0] *= 128.*phi0*xRef[0]*xRef[1];
277 jac[0] = {xRef[1],xRef[0]};
278 jac[0] *= 64.*phi0*phi0;
279 jac[0] += jac0[0];
280 functor( 4, jac );
281 }
282
283 template< class Point, class Functor >
284 void hessianEach ( const Point &x, Functor functor ) const
285 {
286 DUNE_THROW( NotImplemented, "NotImplemented" );
287 DomainType xRef = coordinate( x );
289 functor( 0, hes );
290 functor( 1, hes );
291 functor( 2, hes );
292 functor( 3, hes );
293 }
294
295 int order () const { return 4; }
296
297 std::size_t size () const { return 5; }
298 };
299
300 template< class FunctionSpace, class GridPart, class Storage >
301 class BubbleElementSpace;
302
303 // BubbleElementSpaceTraits
304 // ----------------------
305
306 template< class FunctionSpace, class GridPart, class Storage >
308 {
310
312 typedef GridPart GridPartType;
313
314 static const int codimension = 0;
315
316 private:
317 typedef typename GridPartType::template Codim< codimension >::EntityType EntityType;
318 public:
319 typedef typename FunctionSpaceType :: ScalarFunctionSpaceType ScalarFunctionSpaceType;
320
321 // defined in shapefunctionset.hh
322 // First check if we have a grid with only one element type
323 // (hybrid // grids are not supported with this space yet)
324 // If it only has one get the topologyId of that element type
325 // (note simplex=0 and cube=2^dim-1)
326 static_assert( Dune::Capabilities::hasSingleGeometryType<typename GridPartType::GridType>::v,
327 "bubble elements only implemented for grids with a single geometry type" );
328 static const unsigned int topologyId =
329 Dune::Capabilities::hasSingleGeometryType<typename GridPartType::GridType>::topologyId;
330 // now choose either the simplex or the cube implementation of the
331 // bubble space
334 typedef typename std::conditional< topologyId == 0,
337 // now extend it to a vector valued funcion space
339 // and add some default functionality
341 // finished with the shape function set
342
343 static const int localBlockSize = FunctionSpaceType::dimRange;
344 static const int polynomialOrder = ScalarShapeFunctionSetType::polynomialOrder;
345
348
349 template <class DiscreteFunction, class Operation = DFCommunicationOperation::Add >
351 {
352 typedef Operation OperationType;
354 };
355 };
356
357 // BubbleElementSpace
358 // ----------------
359
361 template< class FunctionSpace, class GridPart, class Storage = CachingStorage >
363 : public DiscreteFunctionSpaceDefault< BubbleElementSpaceTraits< FunctionSpace, GridPart, Storage > >
365 {
368
369 public:
373
375
379
381
382 // type of local interpolation
384
385 // static const InterfaceType defaultInterface = InteriorBorder_InteriorBorder_Interface;
386 static const InterfaceType defaultInterface = GridPartType::indexSetInterfaceType;
387 static const CommunicationDirection defaultDirection = ForwardCommunication;
388
390 const InterfaceType commInterface = defaultInterface,
391 const CommunicationDirection commDirection = defaultDirection )
392 : BaseType( gridPart, commInterface, commDirection ),
394 {
395 }
396
398 {
399 }
400
401 bool contains ( const int codim ) const
402 {
403 // forward to mapper since this information is held there
404 return blockMapper().contains( codim );
405 }
406
407 bool continuous () const { return true; }
408
410 bool continuous ( const IntersectionType & intersection ) const
411 {
412 // forward to the subsapces
413 return true;
414 }
415
417 int order () const
418 {
419 return polynomialOrder;
420 }
421
423 template<class Entity>
424 int order ( const Entity &entity ) const
425 {
426 return polynomialOrder;
427 }
428
430 template< class EntityType >
432 {
433 return BasisFunctionSetType( entity, ShapeFunctionSetType( entity.geometry().type() ) );
434 }
435
440
441 // Non-interface methods
442
446 {
447 return InterpolationType();
448 }
449
455 [[deprecated]]
457 {
458 return interpolation();
459 }
460
461 protected:
463 };
464
465 template< class FunctionSpace,
466 class GridPart,
467 class Storage,
468 class NewFunctionSpace >
469 struct DifferentDiscreteFunctionSpace< BubbleElementSpace < FunctionSpace, GridPart, Storage >, NewFunctionSpace >
470 {
472 };
473
474 template< class FunctionSpace, class GridPart, class Storage >
476 : public EmptyLocalRestrictProlong< BubbleElementSpace< FunctionSpace, GridPart, Storage > >
477 {
479 public:
481 : BaseType()
482 {}
483 };
484
485 } // namespace Fem
486
487} // namespace Dune
488
489#endif // SPACE_P1BUBBLE_HH
Dune::Fem::Double pow(const Dune::Fem::Double &a, const Dune::Fem::Double &b)
Definition: double.hh:947
Definition: bindguard.hh:11
double pow(const Double &a, const double b)
Definition: double.hh:881
DofMapperCode compile(const RefElement &refElement, const LocalCoefficients &localCoefficients)
Definition: compile.hh:50
static const Point & coordinate(const Point &x)
Definition: coordinate.hh:14
Definition: explicitfieldvector.hh:75
Definition: hybrid.hh:86
interface for local functions
Definition: localfunction.hh:77
Definition: space/basisfunctionset/default.hh:52
Default communication handler for discrete functions.
Definition: defaultcommhandler.hh:29
Definition: discretefunctionspace.hh:133
Traits::BasisFunctionSetType BasisFunctionSetType
type of basis function set of this space
Definition: discretefunctionspace.hh:200
GridPartType::IntersectionType IntersectionType
type of the intersections
Definition: discretefunctionspace.hh:225
This is the class with default implementations for discrete function. The methods not marked with hav...
Definition: discretefunctionspace.hh:628
BaseType::EntityType EntityType
Definition: discretefunctionspace.hh:644
BaseType::BlockMapperType BlockMapperType
Definition: discretefunctionspace.hh:660
BaseType::GridPartType GridPartType
Definition: discretefunctionspace.hh:640
GridPartType & gridPart() const
Definition: discretefunctionspace.hh:745
A vector valued function space.
Definition: functionspace.hh:60
FunctionSpaceTraits::RangeType RangeType
Type of range vector (using type of range field) has a Dune::FieldVector type interface.
Definition: functionspaceinterface.hh:71
@ dimDomain
dimension of domain vector space
Definition: functionspaceinterface.hh:46
@ dimRange
dimension of range vector space
Definition: functionspaceinterface.hh:48
FunctionSpaceTraits::DomainType DomainType
Type of domain vector (using type of domain field) has a Dune::FieldVector type interface.
Definition: functionspaceinterface.hh:67
Definition: common/localrestrictprolong.hh:16
Definition: common/localrestrictprolong.hh:107
Definition: code.hh:18
Definition: indexsetdofmapper.hh:719
Definition: localkey.hh:21
FunctionSpace::RangeType RangeType
Definition: p1bubble.hh:33
LocalBubbleElementInterpolation()
Definition: p1bubble.hh:37
LocalBubbleElementInterpolation(LocalBubbleElementInterpolation &&)=default
void operator()(const LocalFunction &lf, LocalDofVector &ldv) const
[Evaluation of local interpolations]
Definition: p1bubble.hh:51
void unbind()
Definition: p1bubble.hh:66
static const int dimRange
Definition: p1bubble.hh:35
void bind(const Entity &)
Definition: p1bubble.hh:65
LocalBubbleElementInterpolation(const LocalBubbleElementInterpolation &)=default
static const int dimDomain
Definition: p1bubble.hh:34
FunctionSpace::DomainType DomainType
Definition: p1bubble.hh:32
Definition: p1bubble.hh:75
std::size_t size() const
[Constructor of LocalKey tripple]
Definition: p1bubble.hh:85
BubbleElementLocalKeyMap(int vertices)
[Constructor of LocalKey tripple]
Definition: p1bubble.hh:77
LocalKey & localKey(std::size_t i)
Definition: p1bubble.hh:87
const LocalKey & localKey(std::size_t i) const
Definition: p1bubble.hh:88
FunctionSpace::DomainType DomainType
Definition: p1bubble.hh:126
void hessianEach(const Point &x, Functor functor) const
Definition: p1bubble.hh:187
SimplexBubbleElementShapeFunctionSet()
Definition: p1bubble.hh:132
FunctionSpace FunctionSpaceType
Definition: p1bubble.hh:125
std::size_t size() const
Definition: p1bubble.hh:200
int order() const
Definition: p1bubble.hh:198
FunctionSpace::JacobianRangeType JacobianRangeType
Definition: p1bubble.hh:129
FunctionSpace::HessianRangeType HessianRangeType
Definition: p1bubble.hh:130
void jacobianEach(const Point &x, Functor functor) const
Definition: p1bubble.hh:159
FunctionSpace::RangeType RangeType
Definition: p1bubble.hh:127
void evaluateEach(const Point &x, Functor functor) const
Definition: p1bubble.hh:142
SimplexBubbleElementShapeFunctionSet(const GeometryType &gt)
Definition: p1bubble.hh:135
void evaluateEach(const Point &x, Functor functor) const
[Main methods for shape functions]
Definition: p1bubble.hh:234
FunctionSpace::HessianRangeType HessianRangeType
Definition: p1bubble.hh:219
Cube2DBubbleElementShapeFunctionSet()
[Main methods for shape functions]
Definition: p1bubble.hh:222
FunctionSpace::RangeType RangeType
Definition: p1bubble.hh:216
int order() const
Definition: p1bubble.hh:295
FunctionSpace::JacobianRangeType JacobianRangeType
Definition: p1bubble.hh:218
void jacobianEach(const Point &x, Functor functor) const
Definition: p1bubble.hh:258
std::size_t size() const
Definition: p1bubble.hh:297
FunctionSpace FunctionSpaceType
Definition: p1bubble.hh:214
FunctionSpace::DomainType DomainType
Definition: p1bubble.hh:215
Cube2DBubbleElementShapeFunctionSet(const GeometryType &gt)
Definition: p1bubble.hh:225
void hessianEach(const Point &x, Functor functor) const
Definition: p1bubble.hh:284
[Class definition for new space]
Definition: p1bubble.hh:365
BaseType::IntersectionType IntersectionType
Definition: p1bubble.hh:377
InterpolationType interpolation() const
return local interpolation object for LocalInterpolation
Definition: p1bubble.hh:445
BlockMapperType & blockMapper() const
obtain the DoF block mapper of this space
Definition: p1bubble.hh:439
BubbleElementSpaceTraits< FunctionSpace, GridPart, Storage > Traits
Definition: p1bubble.hh:370
bool continuous(const IntersectionType &intersection) const
returns true if the space contains only globally continuous functions
Definition: p1bubble.hh:410
int order() const
get global order of space
Definition: p1bubble.hh:417
BubbleElementSpace(GridPartType &gridPart, const InterfaceType commInterface=defaultInterface, const CommunicationDirection commDirection=defaultDirection)
Definition: p1bubble.hh:389
BasisFunctionSetType basisFunctionSet(const EntityType &entity) const
Definition: p1bubble.hh:431
BaseType::BasisFunctionSetType BasisFunctionSetType
Definition: p1bubble.hh:376
static const int polynomialOrder
Definition: p1bubble.hh:372
LocalBubbleElementInterpolation< FunctionSpace > InterpolationType
Definition: p1bubble.hh:383
Traits::ShapeFunctionSetType ShapeFunctionSetType
Definition: p1bubble.hh:371
~BubbleElementSpace()
Definition: p1bubble.hh:397
static const CommunicationDirection defaultDirection
Definition: p1bubble.hh:387
BlockMapperType blockMapper_
Definition: p1bubble.hh:462
bool continuous() const
Definition: p1bubble.hh:407
BaseType::EntityType EntityType
Definition: p1bubble.hh:378
int order(const Entity &entity) const
get global order of space
Definition: p1bubble.hh:424
BaseType::GridPartType GridPartType
Definition: p1bubble.hh:374
InterpolationType interpolation(const EntityType &entity) const
return local interpolation for given entity
Definition: p1bubble.hh:456
BaseType::BlockMapperType BlockMapperType
Definition: p1bubble.hh:380
bool contains(const int codim) const
Definition: p1bubble.hh:401
static const InterfaceType defaultInterface
Definition: p1bubble.hh:386
Definition: p1bubble.hh:308
DefaultBasisFunctionSet< EntityType, ShapeFunctionSetType > BasisFunctionSetType
Definition: p1bubble.hh:340
BubbleElementSpace< FunctionSpace, GridPart, Storage > DiscreteFunctionSpaceType
Definition: p1bubble.hh:309
std::conditional< topologyId==0, ScalarSimplexShapeFunctionSetType, ScalarCubeShapeFunctionSetType >::type ScalarShapeFunctionSetType
Definition: p1bubble.hh:336
FunctionSpace FunctionSpaceType
Definition: p1bubble.hh:311
SimplexBubbleElementShapeFunctionSet< ScalarFunctionSpaceType > ScalarSimplexShapeFunctionSetType
Definition: p1bubble.hh:332
static const int polynomialOrder
Definition: p1bubble.hh:344
FunctionSpaceType::ScalarFunctionSpaceType ScalarFunctionSpaceType
Definition: p1bubble.hh:319
GridPart GridPartType
Definition: p1bubble.hh:312
VectorialShapeFunctionSet< ScalarShapeFunctionSetType, typename FunctionSpaceType::RangeType > ShapeFunctionSetType
Definition: p1bubble.hh:338
IndexSetDofMapper< GridPartType > BlockMapperType
Definition: p1bubble.hh:346
Hybrid::IndexRange< int, FunctionSpaceType::dimRange > LocalBlockIndices
Definition: p1bubble.hh:347
Cube2DBubbleElementShapeFunctionSet< ScalarFunctionSpaceType > ScalarCubeShapeFunctionSetType
Definition: p1bubble.hh:333
Operation OperationType
Definition: p1bubble.hh:352
DefaultCommunicationHandler< DiscreteFunction, Operation > Type
Definition: p1bubble.hh:353
BubbleElementSpace< NewFunctionSpace, GridPart, Storage > Type
Definition: p1bubble.hh:471
DefaultLocalRestrictProlong(const BubbleElementSpace< FunctionSpace, GridPart, Storage > &space)
Definition: p1bubble.hh:480
Definition: shapefunctionset/vectorial.hh:447