3#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEDGBASIS_HH
4#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_LAGRANGEDGBASIS_HH
7#include <dune/common/exceptions.hh>
9#include <dune/functions/functionspacebases/nodes.hh>
10#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
11#include <dune/functions/functionspacebases/flatmultiindex.hh>
12#include <dune/functions/functionspacebases/lagrangebasis.hh>
34template<
typename GV,
int k>
35using LagrangeDGNode = LagrangeNode<GV, k>;
37template<
typename GV,
int k,
class MI>
38class LagrangeDGNodeIndexSet;
41template<
typename GV,
int k,
class MI>
42class LagrangeDGPreBasis
44 static const int dim = GV::dimension;
50 using size_type = std::size_t;
54 const static int dofsPerEdge = k+1;
55 const static int dofsPerTriangle = (k+1)*(k+2)/2;
56 const static int dofsPerQuad = (k+1)*(k+1);
57 const static int dofsPerTetrahedron = (k+1)*(k+2)*(k+3)/6;
58 const static int dofsPerPrism = (k+1)*(k+1)*(k+2)/2;
59 const static int dofsPerHexahedron = (k+1)*(k+1)*(k+1);
60 const static int dofsPerPyramid = (k+1)*(k+2)*(2*k+3)/6;
63 using Node = LagrangeDGNode<GV, k>;
65 using IndexSet = LagrangeDGNodeIndexSet<GV, k, MI>;
68 using MultiIndex = MI;
70 using SizePrefix = Dune::ReservedVector<size_type, 1>;
73 LagrangeDGPreBasis(
const GridView& gv) :
78 void initializeIndices()
88 quadrilateralOffset_ = dofsPerTriangle * gridView_.size(Dune::GeometryTypes::triangle);
93 prismOffset_ = dofsPerTetrahedron * gridView_.size(Dune::GeometryTypes::tetrahedron);
95 hexahedronOffset_ = prismOffset_ + dofsPerPrism * gridView_.size(Dune::GeometryTypes::prism);
97 pyramidOffset_ = hexahedronOffset_ + dofsPerHexahedron * gridView_.size(Dune::GeometryTypes::hexahedron);
105 const GridView& gridView()
const
110 void update(
const GridView& gv)
118 Node makeNode()
const
129 IndexSet makeIndexSet()
const
131 return IndexSet{*
this};
134 size_type size()
const
139 return dofsPerEdge*gridView_.size(0);
142 return dofsPerTriangle*gridView_.size(Dune::GeometryTypes::triangle) + dofsPerQuad*gridView_.size(Dune::GeometryTypes::quadrilateral);
146 return dofsPerTetrahedron*gridView_.size(Dune::GeometryTypes::tetrahedron)
147 + dofsPerPyramid*gridView_.size(Dune::GeometryTypes::pyramid)
148 + dofsPerPrism*gridView_.size(Dune::GeometryTypes::prism)
149 + dofsPerHexahedron*gridView_.size(Dune::GeometryTypes::hexahedron);
152 DUNE_THROW(Dune::NotImplemented,
"No size method for " << dim <<
"d grids available yet!");
156 size_type size(
const SizePrefix prefix)
const
158 assert(prefix.size() == 0 || prefix.size() == 1);
159 return (prefix.size() == 0) ? size() : 0;
163 size_type dimension()
const
168 size_type maxNodeSize()
const
170 return StaticPower<(k+1),GV::dimension>
::power;
176 size_t quadrilateralOffset_;
177 size_t pyramidOffset_;
179 size_t hexahedronOffset_;
184template<
typename GV,
int k,
class MI>
185class LagrangeDGNodeIndexSet
188 static const int dim = GV::dimension;
192 using size_type = std::size_t;
195 using MultiIndex = MI;
197 using PreBasis = LagrangeDGPreBasis<GV, k, MI>;
199 using Node = LagrangeDGNode<GV, k>;
201 LagrangeDGNodeIndexSet(
const PreBasis& preBasis) :
210 void bind(
const Node& node)
224 size_type size()
const
226 return node_->finiteElement().size();
230 template<
typename It>
231 It indices(It it)
const
233 const auto& gridIndexSet = preBasis_->gridView().indexSet();
234 const auto& element = node_->element();
236 for (size_type i = 0, end = size() ; i < end ; ++i, ++it)
242 *it = {preBasis_->dofsPerEdge*gridIndexSet.subIndex(element,0,0) + i};
247 if (element.type().isTriangle())
249 *it = {preBasis_->dofsPerTriangle*gridIndexSet.subIndex(element,0,0) + i};
252 else if (element.type().isQuadrilateral())
254 *it = { preBasis_->quadrilateralOffset_ + preBasis_->dofsPerQuad*gridIndexSet.subIndex(element,0,0) + i};
258 DUNE_THROW(Dune::NotImplemented,
"2d elements have to be triangles or quadrilaterals");
262 if (element.type().isTetrahedron())
264 *it = {preBasis_->dofsPerTetrahedron*gridIndexSet.subIndex(element,0,0) + i};
267 else if (element.type().isPrism())
269 *it = { preBasis_->prismOffset_ + preBasis_->dofsPerPrism*gridIndexSet.subIndex(element,0,0) + i};
272 else if (element.type().isHexahedron())
274 *it = { preBasis_->hexahedronOffset_ + preBasis_->dofsPerHexahedron*gridIndexSet.subIndex(element,0,0) + i};
277 else if (element.type().isPyramid())
279 *it = { preBasis_->pyramidOffset_ + preBasis_->dofsPerPyramid*gridIndexSet.subIndex(element,0,0) + i};
283 DUNE_THROW(Dune::NotImplemented,
"3d elements have to be tetrahedrons, prisms, hexahedrons or pyramids");
286 DUNE_THROW(Dune::NotImplemented,
"No index method for " << dim <<
"d grids available yet!");
292 const PreBasis* preBasis_;
300namespace BasisFactory {
304template<std::
size_t k>
305class LagrangeDGPreBasisFactory
308 static const std::size_t requiredMultiIndexSize = 1;
310 template<
class MultiIndex,
class Gr
idView>
311 auto makePreBasis(
const GridView& gridView)
const
313 return LagrangeDGPreBasis<GridView, k, MultiIndex>(gridView);
329template<std::
size_t k>
332 return Imp::LagrangeDGPreBasisFactory<k>();
350template<
typename GV,
int k>
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:47
auto power(ChildPreBasisFactory &&childPreBasisFactory, const IndexMergingStrategy &ims)
Create a pre-basis factory that can build a PowerPreBasis.
Definition: powerbasis.hh:483
auto lagrangeDG()
Create a pre-basis factory that can create a LagrangeDG pre-basis.
Definition: lagrangedgbasis.hh:330
Definition: polynomial.hh:10