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