dune-fem  2.8-git
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 
23 namespace Dune
24 {
25 
26  namespace Fem
27  {
28 
29  template< class FunctionSpace >
31  {
34  static const int dimDomain = FunctionSpace::dimDomain;
35  static const int dimRange = FunctionSpace::dimRange;
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 
128  static_assert( RangeType::dimension == 1, "This is a scalar shapefunction set" );
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 );
191  HessianRangeType hes;
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 
217  static_assert( RangeType::dimension == 1, "This is a scalar shapefunction set" );
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 );
288  HessianRangeType hes;
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
Definition: bindguard.hh:11
double pow(const Double &a, const double b)
Definition: double.hh:881
static const Point & coordinate(const Point &x)
Definition: coordinate.hh:14
DofMapperCode compile(const RefElement &refElement, const LocalCoefficients &localCoefficients)
Definition: compile.hh:50
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
GridPartType ::IntersectionType IntersectionType
type of the intersections
Definition: discretefunctionspace.hh:225
Traits ::BasisFunctionSetType BasisFunctionSetType
type of basis function set of this space
Definition: discretefunctionspace.hh:200
This is the class with default implementations for discrete function. The methods not marked with hav...
Definition: discretefunctionspace.hh:628
BaseType::BlockMapperType BlockMapperType
Definition: discretefunctionspace.hh:660
BaseType ::EntityType EntityType
Definition: discretefunctionspace.hh:644
BaseType ::GridPartType GridPartType
Definition: discretefunctionspace.hh:640
A vector valued function space.
Definition: functionspace.hh:60
ExplicitFieldVector< FieldMatrix< RangeFieldType, dimDomain, dimDomain >, dimRange > HessianRangeType
Intrinsic type used for the hessian values has a Dune::FieldMatrix type interface.
Definition: functionspaceinterface.hh:79
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::LinearMappingType JacobianRangeType
Intrinsic type used for the jacobian values has a Dune::FieldMatrix type interface.
Definition: functionspaceinterface.hh:75
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:655
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
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 ::HessianRangeType HessianRangeType
Definition: p1bubble.hh:130
FunctionSpace ::DomainType DomainType
Definition: p1bubble.hh:126
FunctionSpace ::RangeType RangeType
Definition: p1bubble.hh:127
void jacobianEach(const Point &x, Functor functor) const
Definition: p1bubble.hh:159
FunctionSpace ::JacobianRangeType JacobianRangeType
Definition: p1bubble.hh:128
void evaluateEach(const Point &x, Functor functor) const
Definition: p1bubble.hh:142
SimplexBubbleElementShapeFunctionSet(const GeometryType &gt)
Definition: p1bubble.hh:135
FunctionSpace ::RangeType RangeType
Definition: p1bubble.hh:216
void evaluateEach(const Point &x, Functor functor) const
[Main methods for shape functions]
Definition: p1bubble.hh:234
Cube2DBubbleElementShapeFunctionSet()
[Main methods for shape functions]
Definition: p1bubble.hh:222
FunctionSpace ::JacobianRangeType JacobianRangeType
Definition: p1bubble.hh:217
int order() const
Definition: p1bubble.hh:295
void jacobianEach(const Point &x, Functor functor) const
Definition: p1bubble.hh:258
std::size_t size() const
Definition: p1bubble.hh:297
FunctionSpace ::HessianRangeType HessianRangeType
Definition: p1bubble.hh:219
FunctionSpace FunctionSpaceType
Definition: p1bubble.hh:214
Cube2DBubbleElementShapeFunctionSet(const GeometryType &gt)
Definition: p1bubble.hh:225
FunctionSpace ::DomainType DomainType
Definition: p1bubble.hh:215
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
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
get basis function set for given entity
Definition: p1bubble.hh:431
BlockMapperType & blockMapper() const
obtain the DoF block mapper of this space
Definition: p1bubble.hh:439
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
std::conditional< topologyId==0, ScalarSimplexShapeFunctionSetType, ScalarCubeShapeFunctionSetType >::type ScalarShapeFunctionSetType
Definition: p1bubble.hh:336
BubbleElementSpace< FunctionSpace, GridPart, Storage > DiscreteFunctionSpaceType
Definition: p1bubble.hh:309
FunctionSpace FunctionSpaceType
Definition: p1bubble.hh:311
SimplexBubbleElementShapeFunctionSet< ScalarFunctionSpaceType > ScalarSimplexShapeFunctionSetType
Definition: p1bubble.hh:332
FunctionSpaceType ::ScalarFunctionSpaceType ScalarFunctionSpaceType
Definition: p1bubble.hh:319
static const int polynomialOrder
Definition: p1bubble.hh:344
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