dune-localfunctions  2.2.1
lobattopoints.hh
Go to the documentation of this file.
1 #ifndef DUNE_LOBATTOBASIS_HH
2 #define DUNE_LOBATTOBASIS_HH
3 
4 #include <fstream>
5 
6 #include <dune/common/forloop.hh>
7 
8 #include <dune/geometry/topologyfactory.hh>
9 #include <dune/geometry/quadraturerules/lobattoquadrature.hh>
10 #include <dune/geometry/referenceelements.hh>
11 #include <dune/geometry/genericgeometry/referenceelements.hh>
12 
16 
17 namespace Dune
18 {
19 
20  template< class Field >
22  {
23  LobattoPoints ( unsigned int order )
24  : points_( 0 )
25  {
26  if( order < 2 )
27  return;
28  points_.resize(order-1);
29 #if HAVE_ALGLIB
30  typedef amp::ampf< Precision< Field >::value > MPField;
31 #else
32  typedef Field MPField;
33 #endif
34  GenericGeometry::LobattoPoints<MPField> lobatto(order+1);
35 
36  for (unsigned int i=1;i<order;++i) {
37  points_[i-1] = field_cast<Field>(lobatto[i].position());
38  }
39  }
40 
41  const unsigned int size()
42  {
43  return points_.size();
44  }
45  const unsigned int order()
46  {
47  return points_.size()+1;
48  }
49  const Field &point(int i)
50  {
51  return points_[i];
52  }
53  std::vector<Field> points_;
54  };
55 
56  template <class Field,class Topology>
58 
59  template <class Field>
60  struct LobattoInnerPoints<Field,GenericGeometry::Point>
61  {
62  static const unsigned int dimension = 0;
63  static unsigned int size(const unsigned int order)
64  {
65  return 1;
66  }
67  template <unsigned int dim>
68  static unsigned int setupSimplex(
69  const unsigned int iup,
70  const unsigned int dimStart,
71  unsigned int startPos,
72  const std::vector<Field> &points1D,
74  {
75  const unsigned int order = points1D.size()+1;
76  unsigned int i = order+1-iup;
77  if (i-2>=points1D.size())
78  {
79  return startPos;
80  }
81  for ( unsigned int d=0;d<dimStart;++d )
82  {
83  points[startPos].point_[d] = -points1D[i-2];
84  }
85 
86  return startPos+1;
87  }
88  template <unsigned int dim>
89  static void setup(const std::vector<Field> &points1D,
91  {
92  points->point_[0] = Zero<Field>();
93  }
94  };
95  template <class Field,class Base>
96  struct LobattoInnerPoints<Field,GenericGeometry::Pyramid<Base> >
97  {
99  static const unsigned int dimension = Base::dimension+1;
100  static unsigned int size(const unsigned int order)
101  {
102  if (order<=dimension) return 0;
103  unsigned int size=0;
104  for (unsigned int o=0;o<order-dimension;++o)
105  size += LobattoBase::size(o+dimension);
106  return size;
107  }
108  template <unsigned int dim>
109  static unsigned int setupSimplex(
110  const unsigned int iup,
111  const unsigned int dimStart,
112  unsigned int startPos,
113  const std::vector<Field> &points1D,
115  {
116  const unsigned int order = points1D.size()+1;
117  unsigned int endPos = startPos;
118  for (unsigned int i=2;i<=order-iup;++i)
119  {
120  endPos = LobattoBase::template setupSimplex<dim>(iup+i-1,dimStart,startPos,points1D,points);
121  for (unsigned int j=startPos;j<endPos;++j)
122  {
123  for (unsigned int d=0;d<dimStart;++d)
124  {
125  if ( d==dimStart-dimension )
126  points[j].point_[d] += dimStart*points1D[i-2];
127  else
128  points[j].point_[d] -= points1D[i-2];
129  }
130  }
131  startPos = endPos;
132  }
133  return endPos;
134  }
135  template <unsigned int dim>
136  static void setup(const std::vector<Field> &points1D,
138  {
139  const unsigned int order = points1D.size()+1;
140  unsigned int startPos=0,endPos=0;
141  for (unsigned int i=2;i<=order;++i)
142  {
143  endPos = LobattoBase::template setupSimplex<dim>(i-1,dimension,startPos,points1D,points);
144 
145  for (unsigned int j=startPos;j<endPos;++j)
146  {
147  for (unsigned int d=0;d<dimension;++d)
148  {
149  if ( d==0 )
150  points[j].point_[d] += 1.+int(dimension)*points1D[i-2];
151  else
152  points[j].point_[d] += 1.-points1D[i-2];
153  points[j].point_[d] /= (dimension+1);
154  }
155  }
156  startPos = endPos;
157  }
158  }
159  };
160  template <class Field,class Base>
161  struct LobattoInnerPoints<Field,GenericGeometry::Prism<Base> >
162  {
164  static const unsigned int dimension = Base::dimension+1;
165  static unsigned int size(const unsigned int order)
166  {
167  return LobattoBase::size(order)*(order-1);
168  }
169  template <unsigned int dim>
170  static unsigned int setupSimplex(
171  const unsigned int iup,
172  const unsigned int dimStart,
173  unsigned int startPos,
174  const std::vector<Field> &points1D,
176  {
177  const unsigned int order = points1D.size()+1;
178  unsigned int endPos = startPos;
179  for (unsigned int i=2;i<=order-iup;++i)
180  {
181  endPos = LobattoBase::template setupSimplex<dim>(iup+i-1,dimStart,startPos,points1D,points);
182  for (unsigned int j=startPos;j<endPos;++j)
183  {
184  for (unsigned int d=0;d<dimStart;++d)
185  {
186  if ( d==dimStart-dimension )
187  points[j].point_[d] += dimStart*points1D[i-2];
188  else
189  points[j].point_[d] -= points1D[i-2];
190  }
191  }
192  startPos = endPos;
193  }
194  return endPos;
195  }
196  template <unsigned int dim>
197  static void setup(const std::vector<Field> &points1D,
199  {
200  const unsigned int order = points1D.size()+1;
201  assert(dim>=dimension);
202  assert(points1D.size()==order-1);
203  LobattoBase::template setup<dim>(points1D,points);
204  const unsigned int baseSize = LobattoBase::size(order);
205  for (unsigned int q=0;q<points1D.size();q++)
206  {
207  for (unsigned int i=0;i<baseSize;++i)
208  {
209  const unsigned int pos = q*baseSize+i;
210  for (unsigned int d=0;d<dimension-1;++d)
211  points[pos].point_[d] = points[i].point_[d];
212  points[pos].point_[dimension-1]=points1D[q];
213  }
214  }
215  }
216  };
217 
218  template< class F, unsigned int dim >
219  struct LobattoPointSet : public EmptyPointSet<F,dim>
220  {
221  // friend class LagrangeCoefficientsFactory<LobattoPointSet,dim,F>;
222  static const unsigned int dimension = dim;
223  typedef F Field;
225  typedef typename Base::LagrangePoint Point;
226  LobattoPointSet(unsigned int order)
227  : Base(order)
228  {
229  }
230 
231  template< class Topology >
232  bool build ( )
233  {
234  unsigned int order = Base::order();
235  LobattoPoints<Field> points1D(order);
236  ForLoop<Setup<Topology>::template InitCodim,0,dimension>::
237  apply(order,points1D.points_,points_);
238  return true;
239  }
240  template< class Topology >
241  static bool supports ( unsigned int order )
242  {
245  return true;
246  else
247  return false;
248  }
249  protected:
250  using Base::points_;
251  private:
252  template <class Topology>
253  struct Setup
254  {
255  template <int pdim>
256  struct InitCodim
257  {
258  static const unsigned int codim = dimension-pdim;
259  static void apply(const unsigned int order,
260  const std::vector<Field> &points1D,
261  std::vector<Point> &points)
262  {
264  ForLoop<InitSub,0,size-1>::apply(order,points1D,points);
265  }
266  template <int i>
267  struct InitSub
268  {
269  typedef typename GenericGeometry::SubTopology<Topology,codim,i>::type SubTopology;
270  static void apply(const unsigned int order,
271  const std::vector<Field> &points1D,
272  std::vector<Point> &points)
273  {
274  Setup<Topology>::template Init<SubTopology>::template apply<i>(order,points1D,points);
275  }
276  };
277  };
278  template <class SubTopology>
279  struct Init
280  {
281  static const unsigned int codimension = dimension - SubTopology::dimension;
282 
283  typedef GenericReferenceElements< Field, dimension > RefElements;
284  typedef GenericReferenceElement< Field, dimension > RefElement;
285  typedef typename RefElement::template Codim< codimension >::Mapping Mapping;
286 
288  template <unsigned int subEntity>
289  static void apply(const unsigned int order,
290  const std::vector<Field> &points1D,
291  std::vector<Point> &points)
292  {
293  unsigned int oldSize = points.size();
294  unsigned int size = InnerPoints::size(order);
295  if (size==0)
296  return;
297  points.resize(oldSize+size);
298  std::vector< LagrangePoint<Field,dimension-codimension> > subPoints(size);
299 
300  InnerPoints::template setup<dimension-codimension>( points1D,&(subPoints[0]) );
301 
302  const GeometryType geoType( Topology::id, dimension );
303  const RefElement &refElement = RefElements::general( geoType );
304  const Mapping &mapping = refElement.template mapping< codimension >( subEntity );
305 
306  LagrangePoint<Field,dimension> *p = &(points[oldSize]);
307  for ( unsigned int nr = 0; nr<size; ++nr, ++p)
308  {
309  p->point_ = mapping.global( subPoints[nr].point_ );
310  p->localKey_ = LocalKey( subEntity, codimension, nr );
311  #ifndef NDEBUG
312  bool test = GenericGeometry::ReferenceElement<Topology,Field>::checkInside(p->point_);
313  if (!test)
314  std::cout << "not inside" << std::endl;
315  #endif
316  }
317  }
318  };
319  };
320 
321  };
322 }
323 #endif // DUNE_LOBATTOBASIS_HH
324