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