dune-localfunctions  2.2.1
equidistantpoints.hh
Go to the documentation of this file.
1 #ifndef DUNE_EQUIDISTANTPOINTS_HH
2 #define DUNE_EQUIDISTANTPOINTS_HH
3 
4 #include <vector>
5 
6 #include <dune/common/fvector.hh>
7 #include <dune/common/forloop.hh>
8 #include <dune/geometry/topologyfactory.hh>
9 #include <dune/geometry/genericgeometry/topologytypes.hh>
10 #include <dune/geometry/genericgeometry/subtopologies.hh>
11 
14 
15 namespace Dune
16 {
17 
18  // Internal Forward Declarations
19  // -----------------------------
20 
21  template< class F, unsigned int dim >
23 
24  // EquidistantPointSetImpl
25  // ----------------------------
26 
27  template< class Topology, class F >
29 
30  template< class F >
31  class EquidistantPointSetImpl< GenericGeometry::Point, F >
32  {
34 
35  typedef GenericGeometry::Point Topology;
36 
37  friend class EquidistantPointSet< F, Topology::dimension >;
38  friend class EquidistantPointSetImpl< GenericGeometry::Prism< Topology >, F >;
39  friend class EquidistantPointSetImpl< GenericGeometry::Pyramid< Topology >, F >;
40 
41  public:
42  typedef F Field;
43 
44  static const unsigned int dimension = Topology::dimension;
45 
46  static unsigned int size ( const unsigned int order )
47  {
48  return 1;
49  }
50 
51  private:
52  template< unsigned int codim, unsigned int dim >
53  static unsigned int setup ( const unsigned int order,
54  unsigned int *count,
56  {
57  assert( codim == 0 );
58  points->localKey_ = LocalKey( 0, 0, count[ 0 ]++ );
59  points->point_ = Field( 0 );
60  return 1;
61  }
62  };
63 
64  template< class BaseTopology, class F >
65  class EquidistantPointSetImpl< GenericGeometry::Prism< BaseTopology >, F >
66  {
68 
69  typedef GenericGeometry::Prism< BaseTopology > Topology;
70 
71  friend class EquidistantPointSet< F, Topology::dimension >;
72  friend class EquidistantPointSetImpl< GenericGeometry::Prism< Topology >, F >;
73  friend class EquidistantPointSetImpl< GenericGeometry::Pyramid< Topology >, F >;
74 
75  typedef EquidistantPointSetImpl< BaseTopology, F > BaseImpl;
76 
77  public:
78  typedef F Field;
79 
80  static const unsigned int dimension = Topology::dimension;
81 
82  static unsigned int size ( const unsigned int order )
83  {
84  return BaseImpl::size( order ) * (order+1);
85  }
86 
87  // private:
88  template< unsigned int codim, unsigned int dim >
89  static unsigned int setup ( const unsigned int order,
90  unsigned int *count,
92  {
93  unsigned int size = 0;
94  unsigned int numBaseN = 0;
95 
96  if( codim < dimension )
97  {
98  const unsigned int vcodim = (codim < dimension ? codim : dimension-1);
100  for( unsigned int i = 1; i < order; ++i )
101  {
102  const unsigned int n = BaseImpl::template setup< vcodim, dim >( order, count, points );
103  for( unsigned int j = 0; j < n; ++j )
104  {
105  LocalKey &key = points->localKey_;
106  key = LocalKey( key.subEntity(), codim, key.index() );
107  points->point_[ dimension-1 ] = Field( i ) / Field( order );
108  ++points;
109  }
110  size += n;
111  }
112  }
113 
114  if( codim > 0 )
115  {
116  const unsigned int vcodim = (codim > 0 ? codim : 1);
117  const unsigned int numBaseM = GenericGeometry::Size< BaseTopology, vcodim-1 >::value;
118  const unsigned int n = BaseImpl::template setup< vcodim-1, dim >( order, count+numBaseN, points );
119  for( unsigned int j = 0; j < n; ++j )
120  {
121  LocalKey &key = points[ j ].localKey_;
122  key = LocalKey( key.subEntity() + numBaseN, codim, key.index() );
123  points[ j + n ].point_ = points[ j ].point_;
124  points[ j + n ].point_[ dimension-1 ] = Field( 1 );
125  points[ j + n ].localKey_ = LocalKey( key.subEntity() + numBaseM, codim, key.index() );
126  ++count[ key.subEntity() + numBaseM ];
127  }
128  size += 2*n;
129  }
130  return size;
131  }
132  };
133 
134  template< class BaseTopology, class F >
135  class EquidistantPointSetImpl< GenericGeometry::Pyramid< BaseTopology >, F >
136  {
138 
139  typedef GenericGeometry::Pyramid< BaseTopology > Topology;
140 
141  friend class EquidistantPointSet< F, Topology::dimension >;
142  friend class EquidistantPointSetImpl< GenericGeometry::Prism< Topology >, F >;
143  friend class EquidistantPointSetImpl< GenericGeometry::Pyramid< Topology >, F >;
144 
145  typedef EquidistantPointSetImpl< BaseTopology, F > BaseImpl;
146 
147  public:
148  typedef F Field;
149 
150  static const unsigned int dimension = Topology::dimension;
151 
152  static unsigned int size ( const unsigned int order )
153  {
154  unsigned int size = BaseImpl::size( order );
155  for( unsigned int i = 1; i <= order; ++i )
156  size += BaseImpl::size( order - i );
157  return size;
158  }
159 
160  // private:
161  template< unsigned int codim, unsigned int dim >
162  static unsigned int setup ( const unsigned int order,
163  unsigned int *count,
165  {
166  unsigned int size = 0;
167  unsigned int numBaseM = 0;
168 
169  if( codim > 0 )
170  {
171  const unsigned int vcodim = (codim > 0 ? codim : 1);
172  numBaseM = GenericGeometry::Size< BaseTopology, vcodim-1 >::value;
173  size = BaseImpl::template setup< vcodim-1, dim >( order, count, points );
174  LagrangePoint< Field, dim > *const end = points + size;
175  for( ; points != end; ++points )
176  {
177  LocalKey &key = points->localKey_;
178  key = LocalKey( key.subEntity(), codim, key.index() );
179  }
180  }
181 
182  if( codim < dimension )
183  {
184  const unsigned int vcodim = (codim < dimension ? codim : dimension-1);
185  for( unsigned int i = order-1; i > 0; --i )
186  {
187  const unsigned int n = BaseImpl::template setup< vcodim, dim >( i, count+numBaseM, points );
188  LagrangePoint< Field, dim > *const end = points + n;
189  for( ; points != end; ++points )
190  {
191  LocalKey &key = points->localKey_;
192  key = LocalKey( key.subEntity()+numBaseM, codim, key.index() );
193  for( unsigned int j = 0; j < dimension-1; ++j )
194  points->point_[ j ] *= Field( i ) / Field( order );
195  points->point_[ dimension-1 ] = Field( order - i ) / Field( order );
196  }
197  size += n;
198  }
199  }
200  else
201  {
202  points->localKey_ = LocalKey( numBaseM, dimension, count[ numBaseM ]++ );
203  points->point_ = Field( 0 );
204  points->point_[ dimension-1 ] = 1;
205  ++size;
206  }
207 
208  return size;
209  }
210  };
211 
212 
213 
214  // LagnragePoints
215  // --------------
216 
217  template< class F, unsigned int dim >
218  class EquidistantPointSet : public EmptyPointSet<F,dim>
219  {
220  template< class T >
221  struct Topology;
222  typedef EmptyPointSet<F,dim> Base;
223  public:
224  static const unsigned int dimension = dim;
225 
226  EquidistantPointSet( unsigned int order )
227  : Base(order)
228  {
229  }
230 
231  template< class T >
232  bool build ( );
233 
234  template< class T >
235  static bool supports ( unsigned int order )
236  {
237  return true;
238  }
239  private:
240  using Base::points_;
241  };
242 
243  template< class F, unsigned int dim >
244  template< class T >
245  struct EquidistantPointSet< F, dim >::Topology
246  {
247  typedef EquidistantPointSetImpl< T, F > Impl;
249 
250  template< int pdim >
251  struct Init
252  {
253  static void apply ( const unsigned int order, LagrangePoint *&p )
254  {
255  const unsigned int size = GenericGeometry::Size< T, dimension-pdim >::value;
256  unsigned int count[ size ];
257  for( unsigned int i = 0; i < size; ++i )
258  count[ i ] = 0;
259  p += Impl::template setup< dimension-pdim, dimension >( order, count, p );
260  }
261  };
262  };
263 
264  template< class F, unsigned int dim >
265  template< class T >
267  {
268  unsigned int order = Base::order();
270  typedef typename Topology< T >::Impl Impl;
271  points_.resize( Impl::size( order ) );
272  LagrangePoint *p = &(points_[ 0 ]);
273  ForLoop< Topology< T >::template Init, 0, dimension >::apply( order, p );
274  return true;
275  }
276 }
277 
278 #endif