dune-fem 2.8.0
Loading...
Searching...
No Matches
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
20namespace 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
Dune::Fem::Double abs(const Dune::Fem::Double &a)
Definition: double.hh:942
Dune::Fem::Double pow(const Dune::Fem::Double &a, const Dune::Fem::Double &b)
Definition: double.hh:947
Definition: bindguard.hh:11
static DUNE_PRIVATE void apply(Args &&... args)
Definition: forloop.hh:23