dune-localfunctions  2.8.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/referenceelements.hh>
10 #include <dune/geometry/type.hh>
11 
14 
15 namespace Dune
16 {
17 
18  // numLagrangePoints
19  // -----------------
20 
21  inline std::size_t numLagrangePoints ( const GeometryType& gt, std::size_t order )
22  {
23  const int dim = gt.dim();
24  if( dim > 0 )
25  {
26  const GeometryType baseGeometryType = Impl::getBase( gt );
27  if( gt.isConical() )
28  {
29  std::size_t size = 0;
30  for( unsigned int o = 0; o <= order; ++o )
31  size += numLagrangePoints( baseGeometryType, o );
32  return size;
33  }
34  else
35  return numLagrangePoints( baseGeometryType, order ) * (order+1);
36  }
37  else
38  return 1;
39  }
40 
41  [[deprecated("Use numLagrangePoints(const GeometryType& gt, std::size_t order ) instead.")]]
42  inline std::size_t numLagrangePoints ( unsigned int topologyId, unsigned int dim, std::size_t order )
43  {
44  return numLagrangePoints ( GeometryType(topologyId, dim), order);
45  }
46 
47 
48 
49  // equidistantLagrangePoints
50  // -------------------------
51 
52  template< class ct, unsigned int cdim >
53  inline static unsigned int equidistantLagrangePoints ( const GeometryType& gt, unsigned int codim, std::size_t order, unsigned int *count, LagrangePoint< ct, cdim > *points )
54  {
55  const unsigned int dim = gt.dim();
56  assert( (0 <= codim) && (codim <= dim) && (dim <= cdim) );
57 
58  if( dim > 0 )
59  {
60  const GeometryType baseGeometryType = Impl::getBase( gt );
61  const unsigned int numBaseN = (codim < dim ? Geo::Impl::size( baseGeometryType.id(), baseGeometryType.dim(), codim ) : 0);
62  const unsigned int numBaseM = (codim > 0 ? Geo::Impl::size( baseGeometryType.id(), baseGeometryType.dim(), codim-1 ) : 0);
63 
64  if( gt.isPrismatic() )
65  {
66  unsigned int size = 0;
67  if( codim < dim )
68  {
69  for( unsigned int i = 1; i < order; ++i )
70  {
71  const unsigned int n = equidistantLagrangePoints( baseGeometryType, codim, order, count, points );
72  for( unsigned int j = 0; j < n; ++j )
73  {
74  LocalKey &key = points->localKey_;
75  key = LocalKey( key.subEntity(), codim, key.index() );
76  points->point_[ dim-1 ] = ct( i ) / ct( order );
77  ++points;
78  }
79  size += n;
80  }
81  }
82 
83  if( codim > 0 )
84  {
85  const unsigned int n = equidistantLagrangePoints( baseGeometryType, codim-1, order, count+numBaseN, points );
86  for( unsigned int j = 0; j < n; ++j )
87  {
88  LocalKey &key = points[ j ].localKey_;
89  key = LocalKey( key.subEntity() + numBaseN, codim, key.index() );
90 
91  points[ j + n ].point_ = points[ j ].point_;
92  points[ j + n ].point_[ dim-1 ] = ct( 1 );
93  points[ j + n ].localKey_ = LocalKey( key.subEntity() + numBaseM, codim, key.index() );
94  ++count[ key.subEntity() + numBaseM ];
95  }
96  size += 2*n;
97  }
98 
99  return size;
100  }
101  else
102  {
103  unsigned int size = (codim > 0 ? equidistantLagrangePoints( baseGeometryType, codim-1, order, count, points ) : 0);
104  LagrangePoint< ct, cdim > *const end = points + size;
105  for( ; points != end; ++points )
106  points->localKey_ = LocalKey( points->localKey_.subEntity(), codim, points->localKey_.index() );
107 
108  if( codim < dim )
109  {
110  for( unsigned int i = order-1; i > 0; --i )
111  {
112  const unsigned int n = equidistantLagrangePoints( baseGeometryType, codim, i, count+numBaseM, points );
113  LagrangePoint< ct, cdim > *const end = points + n;
114  for( ; points != end; ++points )
115  {
116  points->localKey_ = LocalKey( points->localKey_.subEntity()+numBaseM, codim, points->localKey_.index() );
117  for( unsigned int j = 0; j < dim-1; ++j )
118  points->point_[ j ] *= ct( i ) / ct( order );
119  points->point_[ dim-1 ] = ct( order - i ) / ct( order );
120  }
121  size += n;
122  }
123  }
124  else
125  {
126  points->localKey_ = LocalKey( numBaseM, dim, count[ numBaseM ]++ );
127  points->point_ = 0;
128  points->point_[ dim-1 ] = ct( 1 );
129  ++size;
130  }
131 
132  return size;
133  }
134  }
135  else
136  {
137  points->localKey_ = LocalKey( 0, 0, count[ 0 ]++ );
138  points->point_ = 0;
139  return 1;
140  }
141  }
142 
143  template< class ct, unsigned int cdim >
144  [[deprecated("Use equidistantLagrangePoints ( GeometryType gt, ... ) instead.")]]
145  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 )
146  {
147  return equidistantLagrangePoints ( GeometryType(topologyId, dim), codim, order, *count, *points );
148  }
149 
150 
151 
152  // EquidistantPointSet
153  // -------------------
154 
155  template< class F, unsigned int dim >
157  : public EmptyPointSet< F, dim >
158  {
160 
161  public:
162  static const unsigned int dimension = dim;
163 
164  using Base::order;
165 
166  EquidistantPointSet ( std::size_t order ) : Base( order ) {}
167 
168  void build ( GeometryType gt )
169  {
170  assert( gt.dim() == dimension );
171  points_.resize( numLagrangePoints( gt, order() ) );
172 
173  typename Base::LagrangePoint *p = points_.data();
174  std::vector< unsigned int > count;
175  for( unsigned int mydim = 0; mydim <= dimension; ++mydim )
176  {
177  count.resize( Geo::Impl::size( gt.id(), dimension, dimension-mydim ) );
178  std::fill( count.begin(), count.end(), 0u );
179  p += equidistantLagrangePoints( gt, dimension-mydim, order(), count.data(), p );
180  }
181  const auto &refElement = referenceElement<F,dimension>(gt);
182  F weight = refElement.volume()/F(double(points_.size()));
183  for (auto &p : points_)
184  p.weight_ = weight;
185  }
186 
187  template< GeometryType::Id geometryId >
188  bool build ()
189  {
190  build( GeometryType( geometryId ) );
191  return true;
192  }
193 
194  bool buildCube ()
195  {
196  return build< GeometryTypes::cube(dim) > ();
197  }
198 
199  static bool supports ( GeometryType gt, std::size_t order ) { return true; }
200  template< GeometryType::Id geometryId>
201  static bool supports ( std::size_t order ) {
202  return supports( GeometryType( geometryId ), order );
203  }
204 
205  private:
206  using Base::points_;
207  };
208 
209 } // namespace Dune
210 
211 #endif // #ifndef DUNE_LOCALFUNCTIONS_LAGRANGE_EQUIDISTANTPOINTS_HH
Definition: bdfmcube.hh:16
std::size_t numLagrangePoints(const GeometryType &gt, std::size_t order)
Definition: equidistantpoints.hh:21
static unsigned int equidistantLagrangePoints(const GeometryType &gt, unsigned int codim, std::size_t order, unsigned int *count, LagrangePoint< ct, cdim > *points)
Definition: equidistantpoints.hh:53
Describe position of one degree of freedom.
Definition: localkey.hh:21
unsigned int index() const
Return offset within subentity.
Definition: localkey.hh:66
unsigned int subEntity() const
Return number of associated subentity.
Definition: localkey.hh:54
Definition: emptypoints.hh:16
Field weight_
Definition: emptypoints.hh:46
Vector point_
Definition: emptypoints.hh:44
LocalKey localKey_
Definition: emptypoints.hh:45
Definition: emptypoints.hh:54
std::size_t order() const
Definition: emptypoints.hh:93
std::vector< LagrangePoint > points_
Definition: emptypoints.hh:105
Definition: equidistantpoints.hh:158
std::size_t order() const
Definition: emptypoints.hh:93
bool build()
Definition: equidistantpoints.hh:188
static bool supports(std::size_t order)
Definition: equidistantpoints.hh:201
static const unsigned int dimension
Definition: equidistantpoints.hh:162
void build(GeometryType gt)
Definition: equidistantpoints.hh:168
static bool supports(GeometryType gt, std::size_t order)
Definition: equidistantpoints.hh:199
bool buildCube()
Definition: equidistantpoints.hh:194
EquidistantPointSet(std::size_t order)
Definition: equidistantpoints.hh:166