1#ifndef DUNE_LOBATTOBASIS_HH
2#define DUNE_LOBATTOBASIS_HH
6#include <dune/geometry/type.hh>
7#include <dune/geometry/topologyfactory.hh>
8#include <dune/geometry/quadraturerules.hh>
9#include <dune/geometry/referenceelements.hh>
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>
24 template <
class Field>
27 template <
class Po
ints1DType>
28 static int size(
const GeometryType gt,
const Points1DType &points1D)
30 if (gt.dim()==0)
return 1;
31 else if (gt.isPrismatic())
32 return Builder<Field>::size(Impl::getBase(gt),points1D)*points1D.size();
35 std::cout <<
"Not implemented for pyramid geometries still missing!\n";
39 template <
unsigned int dim,
class Po
ints1DType>
40 static void setup(
const GeometryType gt,
const Points1DType &points1D,
41 LagrangePoint< Field, dim > *points )
50 points->point_[0] = Zero<Field>();
53 else if (gt.isPrismatic())
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)
61 Field weight = points[i].weight_;
62 for (
unsigned int q=0;q<points1D.size();q++)
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;
74 std::cout <<
"Not implemented for pyramid geometries still missing!\n";
81 template<
class F,
unsigned int dim >
82 struct PointSetFromQuadrature :
public EmptyPointSet<F,dim>
84 static const unsigned int dimension = dim;
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)
93 template <GeometryType::Id geometryId,
class Quad>
94 bool build (
const Quad& quadFactory)
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)
108 Field p = field_cast<Field>(quad[i].position());
112 withEndPoints =
true;
113 vertexWeight = quad[i].weight();
116 points1D.push_back(std::make_pair(p,quad[i].weight()));
120 apply(gt,order,points1D,vertexWeight,points_);
122 Setup::template InitCodim<dimension>::
123 apply(gt,order,points1D,vertexWeight,points_);
126 static bool supports ( GeometryType gt,
int order )
130 template< GeometryType::Id geometryId>
131 static bool supports ( std::size_t order ) {
132 return supports( GeometryType( geometryId ), order );
134 unsigned int quadOrder()
const
140 unsigned int quadOrder_;
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)
153 const unsigned int subEntities = Dune::Geo::Impl::size(gt.id(),gt.dim(),codim);
154 for (
unsigned int subEntity=0;subEntity<subEntities;++subEntity)
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]) );
164 const auto &refElement = referenceElement<Field,dimension>(gt);
165 const auto &mapping = refElement.template geometry< codim >( subEntity );
167 LagrangePoint<Field,dimension> *p = &(points[oldSize]);
168 for (
unsigned int nr = 0; nr<size; ++nr, ++p)
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 );
184 template<
class F,
unsigned int dim >
185 struct GaussLobattoPointSet :
public PointSetFromQuadrature<F,dim>
187 static const unsigned int dimension = dim;
189 typedef PointSetFromQuadrature<F,dim> Base;
190 typedef typename Base::LagrangePoint Point;
193 static const int pointSetId = Dune::QuadratureType::GaussLobatto;
195 GaussLobattoPointSet(
unsigned int order)
198 template< GeometryType::Id geometryId >
202 auto quadFactory = [](
int order)
203 {
return Dune::QuadratureRules<Field,1>::rule(
204 Dune::GeometryTypes::line, pol2QuadOrder(order),
205 Dune::QuadratureType::GaussLobatto);
207 return Base::template build<geometryId>(quadFactory);
209 static unsigned int pol2QuadOrder(
int order)
211 return (order>0)? 2*order-1 : 0;
213 static unsigned int quad2PolOrder(
int order)
218 static auto buildCubeQuadrature(
unsigned int quadOrder)
220 using namespace Impl;
221 GaussLobattoPointSet ps(quad2PolOrder(quadOrder));
222 ps.template build<GeometryTypes::cube(dim)>();
228 template<
class F,
unsigned int dim >
229 struct GaussLegendrePointSet :
public PointSetFromQuadrature<F,dim>
231 static const unsigned int dimension = dim;
233 typedef PointSetFromQuadrature<F,dim> Base;
234 typedef typename Base::LagrangePoint Point;
237 static const int pointSetId = Dune::QuadratureType::GaussLegendre;
239 GaussLegendrePointSet(
unsigned int order)
242 template< GeometryType::Id geometryId >
246 auto quadFactory = [](
int order)
247 {
return Dune::QuadratureRules<Field,1>::rule(
248 Dune::GeometryTypes::line, pol2QuadOrder(order), Dune::QuadratureType::GaussLegendre);
250 return Base::template build<GeometryType(geometryId)>(quadFactory);
253 static unsigned int pol2QuadOrder(
int order)
257 static unsigned int quad2PolOrder(
int order)
262 static auto buildCubeQuadrature(
unsigned int quadOrder)
264 using namespace Impl;
265 GaussLegendrePointSet ps(quad2PolOrder(quadOrder));
266 ps.template build<GeometryTypes::cube(dim)>();
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