dune-fem  2.8-git
quadratureinterpolation.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/geometry/type.hh>
7 #include <dune/geometry/topologyfactory.hh>
8 #include <dune/geometry/quadraturerules.hh>
9 #include <dune/geometry/referenceelements.hh>
10 
12 
13 #if HAVE_DUNE_LOCALFUNCTIONS
14 #include <dune/localfunctions/utility/field.hh>
15 #include <dune/localfunctions/lagrange/lagrangecoefficients.hh>
16 #include <dune/localfunctions/lagrange/emptypoints.hh>
17 #include <dune/localfunctions/lagrange/equidistantpoints.hh>
18 
19 
20 namespace Dune
21 {
22  namespace Impl
23  {
24  template <class Field>
25  struct Builder
26  {
27  template <class Points1DType>
28  static int size(const GeometryType gt, const Points1DType &points1D)
29  {
30  if (gt.dim()==0) return 1;
31  else if (gt.isPrismatic())
32  return Builder<Field>::size(Impl::getBase(gt),points1D)*points1D.size(); // (order-1);
33  else
34  {
35  std::cout << "Not implemented for pyramid geometries still missing!\n";
36  std::abort();
37  }
38  }
39  template <unsigned int dim, class Points1DType>
40  static void setup(const GeometryType gt, const Points1DType &points1D,
41  LagrangePoint< Field, dim > *points )
42  {
43  if (dim==0)
44  {
45  points->weight_ = 1.;
46  return;
47  }
48  if (gt.dim()==0)
49  {
50  points->point_[0] = Zero<Field>();
51  points->weight_ = 1.;
52  }
53  else if (gt.isPrismatic())
54  {
55  GeometryType baseGt = Impl::getBase(gt);
56  assert(dim>=gt.dim());
57  Builder<Field>::template setup<dim>(baseGt,points1D,points);
58  const unsigned int baseSize = Builder::size(baseGt,points1D);
59  for (unsigned int i=0;i<baseSize;++i)
60  {
61  Field weight = points[i].weight_;
62  for (unsigned int q=0;q<points1D.size();q++)
63  {
64  const unsigned int pos = q*baseSize+i;
65  for (unsigned int d=0;d<gt.dim()-1;++d)
66  points[pos].point_[d] = points[i].point_[d];
67  points[pos].point_[gt.dim()-1]=points1D[q].first;
68  points[pos].weight_ = weight*points1D[q].second;
69  }
70  }
71  }
72  else
73  {
74  std::cout << "Not implemented for pyramid geometries still missing!\n";
75  std::abort();
76  }
77  }
78  };
79  } // namespace Impl
80 
81  template< class F, unsigned int dim >
82  struct PointSetFromQuadrature : public EmptyPointSet<F,dim>
83  {
84  static const unsigned int dimension = dim;
85  typedef F Field;
86  typedef EmptyPointSet<F,dim> Base;
87  typedef typename Base::LagrangePoint Point;
88  typedef std::vector<std::pair<Field,Field>> Points1DType;
89  PointSetFromQuadrature(unsigned int order)
90  : Base(order), quadOrder_(-1)
91  {}
92 
93  template <GeometryType::Id geometryId, class Quad>
94  bool build (const Quad& quadFactory)
95  {
96  constexpr GeometryType gt(geometryId);
97  unsigned int order = Base::order();
98  const auto &quad = quadFactory(order);
99  quadOrder_ = quad.order();
100  assert( quad.size() == order+1 );
101  bool withEndPoints = false;
102  Points1DType points1D;
103  Field vertexWeight = 1;
104  for (unsigned int i=0;i<=order;++i)
105  {
106  // remove corner points if part of the quadrature - are added in by
107  // topology construction of points
108  Field p = field_cast<Field>(quad[i].position());
109  Field q = p-1.;
110  if (std::abs(p)<1e-12 || std::abs(q)<1e-12)
111  {
112  withEndPoints = true;
113  vertexWeight = quad[i].weight(); // assuming weight is identical for both end points
114  }
115  else
116  points1D.push_back(std::make_pair(p,quad[i].weight()));
117  }
118  if (withEndPoints)
120  apply(gt,order,points1D,vertexWeight,points_);
121  else
122  Setup::template InitCodim<dimension>::
123  apply(gt,order,points1D,vertexWeight,points_);
124  return true;
125  }
126  static bool supports ( GeometryType gt, int order )
127  {
128  return gt.isCube();
129  }
130  template< GeometryType::Id geometryId>
131  static bool supports ( std::size_t order ) {
132  return supports( GeometryType( geometryId ), order );
133  }
134  unsigned int quadOrder() const
135  {
136  return quadOrder_;
137  }
138  protected:
139  using Base::points_;
140  unsigned int quadOrder_;
141  private:
142  struct Setup
143  {
144  template <int pdim>
145  struct InitCodim
146  {
147  static const unsigned int codim = dimension-pdim;
148  static void apply(GeometryType gt, const unsigned int order,
149  const Points1DType &points1D,
150  const Field &vertexWeight,
151  std::vector<Point> &points)
152  {
153  const unsigned int subEntities = Dune::Geo::Impl::size(gt.id(),gt.dim(),codim);
154  for (unsigned int subEntity=0;subEntity<subEntities;++subEntity)
155  {
156  GeometryType subGt(Impl::baseTopologyId(gt.id(),gt.dim(),codim),gt.dim()-codim);
157  unsigned int oldSize = points.size();
158  unsigned int size = Impl::Builder<Field>::size(subGt,points1D);
159  if (size==0) continue;
160  points.resize(oldSize+size);
161  std::vector< LagrangePoint<Field,dimension-codim> > subPoints(size);
162  Impl::Builder<Field>::template setup<dimension-codim>( subGt, points1D, &(subPoints[0]) );
163 
164  const auto &refElement = referenceElement<Field,dimension>(gt);
165  const auto &mapping = refElement.template geometry< codim >( subEntity );
166 
167  LagrangePoint<Field,dimension> *p = &(points[oldSize]);
168  for ( unsigned int nr = 0; nr<size; ++nr, ++p)
169  {
170  p->point_ = mapping.global( subPoints[nr].point_ );
171  p->weight_ = subPoints[nr].weight_ * std::pow(vertexWeight,codim)*refElement.volume();
172  p->localKey_ = LocalKey( subEntity, codim, nr );
173  }
174  }
175  }
176  };
177  };
178  };
179 
181  // Some pointsets from dune-geometry
183 
184  template< class F, unsigned int dim >
185  struct GaussLobattoPointSet : public PointSetFromQuadrature<F,dim>
186  {
187  static const unsigned int dimension = dim;
188  typedef F Field;
189  typedef PointSetFromQuadrature<F,dim> Base;
190  typedef typename Base::LagrangePoint Point;
191 
192  // enum identifier from dune-geometry QuadratureRules
193  static const int pointSetId = Dune::QuadratureType::GaussLobatto;
194 
195  GaussLobattoPointSet(unsigned int order)
196  : Base(order)
197  {}
198  template< GeometryType::Id geometryId >
199  bool build ()
200  {
201  // get LobattoQuad with order+1 points
202  auto quadFactory = [](int order)
203  { return Dune::QuadratureRules<Field,1>::rule(
204  Dune::GeometryTypes::line, pol2QuadOrder(order),
205  Dune::QuadratureType::GaussLobatto);
206  };
207  return Base::template build<geometryId>(quadFactory);
208  }
209  static unsigned int pol2QuadOrder(int order)
210  {
211  return (order>0)? 2*order-1 : 0;
212  }
213  static unsigned int quad2PolOrder(int order)
214  {
215  return order/2 + 1;
216  }
217 
218  static auto buildCubeQuadrature(unsigned int quadOrder)
219  {
220  using namespace Impl;
221  GaussLobattoPointSet ps(quad2PolOrder(quadOrder));
222  ps.template build<GeometryTypes::cube(dim)>();
223  return ps;
224  }
225  };
226 
227 
228  template< class F, unsigned int dim >
229  struct GaussLegendrePointSet : public PointSetFromQuadrature<F,dim>
230  {
231  static const unsigned int dimension = dim;
232  typedef F Field;
233  typedef PointSetFromQuadrature<F,dim> Base;
234  typedef typename Base::LagrangePoint Point;
235 
236  // enum identifier from dune-geometry QuadratureRules
237  static const int pointSetId = Dune::QuadratureType::GaussLegendre;
238 
239  GaussLegendrePointSet(unsigned int order)
240  : Base(order)
241  {}
242  template< GeometryType::Id geometryId >
243  bool build ()
244  {
245  // get LobattoQuad with order+1 points
246  auto quadFactory = [](int order)
247  { return Dune::QuadratureRules<Field,1>::rule(
248  Dune::GeometryTypes::line, pol2QuadOrder(order), Dune::QuadratureType::GaussLegendre);
249  };
250  return Base::template build<GeometryType(geometryId)>(quadFactory);
251  }
252 
253  static unsigned int pol2QuadOrder(int order)
254  {
255  return 2*order+1;
256  }
257  static unsigned int quad2PolOrder(int order)
258  {
259  return order/2;
260  }
261 
262  static auto buildCubeQuadrature(unsigned int quadOrder)
263  {
264  using namespace Impl;
265  GaussLegendrePointSet ps(quad2PolOrder(quadOrder));
266  ps.template build<GeometryTypes::cube(dim)>();
267  return ps;
268  }
269 
270  };
271 } // namespace DUNE
272 
273 #endif // HAVE_DUNE_LOCALFUNCTIONS
274 
275 #endif // DUNE_LOBATTOBASIS_HH
Definition: bindguard.hh:11
double pow(const Double &a, const double b)
Definition: double.hh:881
Double abs(const Double &a)
Definition: double.hh:871
static DUNE_PRIVATE void apply(Args &&... args)
Definition: forloop.hh:23