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