dune-localfunctions  2.7.0
equidistantpoints.hh
Go to the documentation of this file.
1 #ifndef DUNE_LOCALFUNCTIONS_LAGRANGE_EQUIDISTANTPOINTS_HH
2 #define DUNE_LOCALFUNCTIONS_LAGRANGE_EQUIDISTANTPOINTS_HH
3 
4 #include <cstddef>
5 
6 #include <algorithm>
7 #include <vector>
8 
9 #include <dune/geometry/referenceelementimplementation.hh>
10 #include <dune/geometry/type.hh>
11 
14 
15 namespace Dune
16 {
17 
18  // numLagrangePoints
19  // -----------------
20 
21  inline std::size_t numLagrangePoints ( unsigned int topologyId, int dim, std::size_t order )
22  {
23  assert( topologyId < Impl::numTopologies( dim ) );
24 
25  if( dim > 0 )
26  {
27  const unsigned int baseId = Impl::baseTopologyId( topologyId, dim );
28  if( Impl::isPyramid( topologyId, dim ) )
29  {
30  std::size_t size = 0;
31  for( unsigned int o = 0; o <= order; ++o )
32  size += numLagrangePoints( baseId, dim-1, o );
33  return size;
34  }
35  else
36  return numLagrangePoints( baseId, dim-1, order ) * (order+1);
37  }
38  else
39  return 1;
40  }
41 
42 
43 
44  // equidistantLagrangePoints
45  // -------------------------
46 
47  template< class ct, unsigned int cdim >
48  inline static unsigned int equidistantLagrangePoints ( unsigned int topologyId, unsigned int dim, unsigned int codim, std::size_t order, unsigned int *count, LagrangePoint< ct, cdim > *points )
49  {
50  assert( (0 <= codim) && (codim <= dim) && (dim <= cdim) );
51  assert( topologyId < Impl::numTopologies( dim ) );
52 
53  if( dim > 0 )
54  {
55  const unsigned int baseId = Impl::baseTopologyId( topologyId, dim );
56  const unsigned int numBaseN = (codim < dim ? Geo::Impl::size( baseId, dim-1, codim ) : 0);
57  const unsigned int numBaseM = (codim > 0 ? Geo::Impl::size( baseId, dim-1, codim-1 ) : 0);
58 
59  if( Impl::isPrism( topologyId, dim ) )
60  {
61  unsigned int size = 0;
62  if( codim < dim )
63  {
64  for( unsigned int i = 1; i < order; ++i )
65  {
66  const unsigned int n = equidistantLagrangePoints( baseId, dim-1, codim, order, count, points );
67  for( unsigned int j = 0; j < n; ++j )
68  {
69  LocalKey &key = points->localKey_;
70  key = LocalKey( key.subEntity(), codim, key.index() );
71  points->point_[ dim-1 ] = ct( i ) / ct( order );
72  ++points;
73  }
74  size += n;
75  }
76  }
77 
78  if( codim > 0 )
79  {
80  const unsigned int n = equidistantLagrangePoints( baseId, dim-1, codim-1, order, count+numBaseN, points );
81  for( unsigned int j = 0; j < n; ++j )
82  {
83  LocalKey &key = points[ j ].localKey_;
84  key = LocalKey( key.subEntity() + numBaseN, codim, key.index() );
85 
86  points[ j + n ].point_ = points[ j ].point_;
87  points[ j + n ].point_[ dim-1 ] = ct( 1 );
88  points[ j + n ].localKey_ = LocalKey( key.subEntity() + numBaseM, codim, key.index() );
89  ++count[ key.subEntity() + numBaseM ];
90  }
91  size += 2*n;
92  }
93 
94  return size;
95  }
96  else
97  {
98  unsigned int size = (codim > 0 ? equidistantLagrangePoints( baseId, dim-1, codim-1, order, count, points ) : 0);
99  LagrangePoint< ct, cdim > *const end = points + size;
100  for( ; points != end; ++points )
101  points->localKey_ = LocalKey( points->localKey_.subEntity(), codim, points->localKey_.index() );
102 
103  if( codim < dim )
104  {
105  for( unsigned int i = order-1; i > 0; --i )
106  {
107  const unsigned int n = equidistantLagrangePoints( baseId, dim-1, codim, i, count+numBaseM, points );
108  LagrangePoint< ct, cdim > *const end = points + n;
109  for( ; points != end; ++points )
110  {
111  points->localKey_ = LocalKey( points->localKey_.subEntity()+numBaseM, codim, points->localKey_.index() );
112  for( unsigned int j = 0; j < dim-1; ++j )
113  points->point_[ j ] *= ct( i ) / ct( order );
114  points->point_[ dim-1 ] = ct( order - i ) / ct( order );
115  }
116  size += n;
117  }
118  }
119  else
120  {
121  points->localKey_ = LocalKey( numBaseM, dim, count[ numBaseM ]++ );
122  points->point_ = 0;
123  points->point_[ dim-1 ] = ct( 1 );
124  ++size;
125  }
126 
127  return size;
128  }
129  }
130  else
131  {
132  points->localKey_ = LocalKey( 0, 0, count[ 0 ]++ );
133  points->point_ = 0;
134  return 1;
135  }
136  }
137 
138 
139 
140  // EquidistantPointSet
141  // -------------------
142 
143  template< class F, unsigned int dim >
145  : public EmptyPointSet< F, dim >
146  {
148 
149  public:
150  static const unsigned int dimension = dim;
151 
152  using Base::order;
153 
154  EquidistantPointSet ( std::size_t order ) : Base( order ) {}
155 
156  void build ( GeometryType gt )
157  {
158  assert( gt.dim() == dimension );
159  points_.resize( numLagrangePoints( gt.id(), dimension, order() ) );
160 
161  typename Base::LagrangePoint *p = points_.data();
162  std::vector< unsigned int > count;
163  for( unsigned int mydim = 0; mydim <= dimension; ++mydim )
164  {
165  count.resize( Geo::Impl::size( gt.id(), dimension, dimension-mydim ) );
166  std::fill( count.begin(), count.end(), 0u );
167  p += equidistantLagrangePoints( gt.id(), dimension, dimension-mydim, order(), count.data(), p );
168  }
169  }
170 
171  template< class T >
172  bool build ()
173  {
174  build( GeometryType( T() ) );
175  return true;
176  }
177 
178  template< class T >
179  static bool supports ( std::size_t order ) { return true; }
180 
181  private:
182  using Base::points_;
183  };
184 
185 } // namespace Dune
186 
187 #endif // #ifndef DUNE_LOCALFUNCTIONS_LAGRANGE_EQUIDISTANTPOINTS_HH
Dune::LagrangePoint::point_
Vector point_
Definition: emptypoints.hh:39
Dune::EquidistantPointSet::supports
static bool supports(std::size_t order)
Definition: equidistantpoints.hh:179
Dune::LagrangePoint::localKey_
LocalKey localKey_
Definition: emptypoints.hh:40
Dune
Definition: bdfmcube.hh:15
Dune::EmptyPointSet::order
std::size_t order() const
Definition: emptypoints.hh:87
Dune::EquidistantPointSet::build
bool build()
Definition: equidistantpoints.hh:172
Dune::LocalKey::subEntity
unsigned int subEntity() const
Return number of associated subentity.
Definition: localkey.hh:54
emptypoints.hh
field.hh
Dune::equidistantLagrangePoints
static unsigned int equidistantLagrangePoints(unsigned int topologyId, unsigned int dim, unsigned int codim, std::size_t order, unsigned int *count, LagrangePoint< ct, cdim > *points)
Definition: equidistantpoints.hh:48
Dune::EquidistantPointSet::build
void build(GeometryType gt)
Definition: equidistantpoints.hh:156
Dune::numLagrangePoints
std::size_t numLagrangePoints(unsigned int topologyId, int dim, std::size_t order)
Definition: equidistantpoints.hh:21
Dune::EquidistantPointSet::EquidistantPointSet
EquidistantPointSet(std::size_t order)
Definition: equidistantpoints.hh:154
Dune::LocalKey::index
unsigned int index() const
Return offset within subentity.
Definition: localkey.hh:66
Dune::EmptyPointSet::points_
std::vector< LagrangePoint > points_
Definition: emptypoints.hh:99
Dune::EquidistantPointSet
Definition: equidistantpoints.hh:144
Dune::LagrangePoint
Definition: emptypoints.hh:15
Dune::EquidistantPointSet::dimension
static const unsigned int dimension
Definition: equidistantpoints.hh:150
Dune::EmptyPointSet
Definition: emptypoints.hh:47
Dune::LocalKey
Describe position of one degree of freedom.
Definition: localkey.hh:20