From b6bf445d925f83361ca294646a01915d5866df75 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Wed, 1 May 2019 13:39:28 +0200 Subject: [PATCH 1/3] Added a generic version of the LgrangeBasis that allows to specify the range type of the local basis. --- src/amdis/functions/CMakeLists.txt | 1 + src/amdis/functions/GenericLagrangeBasis.hpp | 511 +++++++++++++++++++ 2 files changed, 512 insertions(+) create mode 100644 src/amdis/functions/GenericLagrangeBasis.hpp diff --git a/src/amdis/functions/CMakeLists.txt b/src/amdis/functions/CMakeLists.txt index 837a5f64..4607552e 100644 --- a/src/amdis/functions/CMakeLists.txt +++ b/src/amdis/functions/CMakeLists.txt @@ -2,6 +2,7 @@ install(FILES FunctionFromCallable.hpp + GenericLagrangeBasis.hpp GlobalIdSet.hpp HierarchicNodeToRangeMap.hpp Interpolate.hpp diff --git a/src/amdis/functions/GenericLagrangeBasis.hpp b/src/amdis/functions/GenericLagrangeBasis.hpp new file mode 100644 index 00000000..9601a84e --- /dev/null +++ b/src/amdis/functions/GenericLagrangeBasis.hpp @@ -0,0 +1,511 @@ +#pragma once + +#include +#include + +#include + +#include +#include +#include + +#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7) +#error "This implementation of GenericLagrange basis requires dune 2.7" +#endif + + +namespace Dune { +namespace Functions { + +/** + * A GenericLagrangeBasis is a GlobalBasis similar to \ref LagrangeBasis, + * that allows to specify the range type of the local basis as template + * parameter. + * + * Construct a GenericLagrangeBasis with \ref lagrange() that now accepts + * two template parameters, k... the polynomial degree and R... the range + * type of the local basis. + * + * Example: + * ``` + * auto basis = makeBasis(gridView, lagrange<1, double>()); + * ``` + **/ + +// forward declarations +template +class GenericLagrangeNode; + +template +class GenericLagrangeNodeIndexSet; + +template +class GenericLagrangePreBasis; + + + +/** + * \brief A pre-basis for a PQ-lagrange bases with given order and range type + * + * \ingroup FunctionSpaceBasesImplementations + * + * \tparam GV The grid view that the FE basis is defined on + * \tparam k The polynomial order of ansatz functions + * \tparam R Range type used for shape function values + * \tparam MI Type to be used for multi-indices + * + * \note This only works for certain grids. The following restrictions hold + * - If k is no larger than 2, then the grids can have any dimension + * - If k is larger than 3 then the grid must be two-dimensional + * - If k is 3, then the grid can be 3d *if* it is a simplex grid + */ +template +class GenericLagrangePreBasis +{ + static const int dim = GV::dimension; + +public: + /// The grid view that the FE basis is defined on + using GridView = GV; + + /// Type used for indices and size information + using size_type = std::size_t; + +private: + template + friend class GenericLagrangeNodeIndexSet; + + // Precompute the number of dofs per entity type + const static size_type dofsPerVertex = + k == 0 ? (dim == 0 ? 1 : 0) : 1; + const static size_type dofsPerEdge = + k == 0 ? (dim == 1 ? 1 : 0) : k-1; + const static size_type dofsPerTriangle = + k == 0 ? (dim == 2 ? 1 : 0) : (k-1)*(k-2)/2; + const static size_type dofsPerQuad = + k == 0 ? (dim == 2 ? 1 : 0) : (k-1)*(k-1); + const static size_type dofsPerTetrahedron = + k == 0 ? (dim == 3 ? 1 : 0) : (k-3)*(k-2)*(k-1)/6; + const static size_type dofsPerPrism = + k == 0 ? (dim == 3 ? 1 : 0) : (k-1)*(k-1)*(k-2)/2; + const static size_type dofsPerHexahedron = + k == 0 ? (dim == 3 ? 1 : 0) : (k-1)*(k-1)*(k-1); + const static size_type dofsPerPyramid = + k == 0 ? (dim == 3 ? 1 : 0) : (k-2)*(k-1)*(2*k-3)/6; + +public: + /// Template mapping root tree path to type of created tree node + using Node = GenericLagrangeNode; + + /// Template mapping root tree path to type of created tree node index set + using IndexSet = GenericLagrangeNodeIndexSet; + + /// Type used for global numbering of the basis vectors + using MultiIndex = MI; + + /// Type used for prefixes handed to the size() method + using SizePrefix = Dune::ReservedVector; + + /// Constructor for a given grid view object + GenericLagrangePreBasis(GridView const& gridView) + : gridView_(gridView) + {} + + /// Initialize the global indices + void initializeIndices() + { + vertexOffset_ = 0; + edgeOffset_ = vertexOffset_ + dofsPerVertex * size_type(gridView_.size(dim)); + + if (dim >= 2) { + triangleOffset_ = edgeOffset_ + dofsPerEdge * size_type(gridView_.size(dim-1)); + quadrilateralOffset_ = triangleOffset_ + dofsPerTriangle * size_type(gridView_.size(Dune::GeometryTypes::triangle)); + } + + if (dim == 3) { + tetrahedronOffset_ = quadrilateralOffset_ + dofsPerQuad * size_type(gridView_.size(Dune::GeometryTypes::quadrilateral)); + prismOffset_ = tetrahedronOffset_ + dofsPerTetrahedron * size_type(gridView_.size(Dune::GeometryTypes::tetrahedron)); + hexahedronOffset_ = prismOffset_ + dofsPerPrism * size_type(gridView_.size(Dune::GeometryTypes::prism)); + pyramidOffset_ = hexahedronOffset_ + dofsPerHexahedron * size_type(gridView_.size(Dune::GeometryTypes::hexahedron)); + } + } + + /// Obtain the grid view that the basis is defined on + GridView const& gridView() const + { + return gridView_; + } + + /// Update the stored grid view, to be called if the grid has changed + void update(GridView const& gridView) + { + gridView_ = gridView; + } + + /// Create tree node + Node makeNode() const + { + return Node{}; + } + + /** + * \brief Create tree node index set + * + * Create an index set suitable for the tree node obtained + * by makeNode(). + */ + IndexSet makeIndexSet() const + { + return IndexSet{*this}; + } + + /// Same as size(prefix) with empty prefix + size_type size() const + { + switch (dim) { + case 1: + { + return dofsPerVertex * size_type(gridView_.size(dim)) + + dofsPerEdge * size_type(gridView_.size(dim-1)); + } + case 2: + { + return dofsPerVertex * size_type(gridView_.size(dim)) + + dofsPerEdge * size_type(gridView_.size(dim-1)) + + dofsPerTriangle * size_type(gridView_.size(Dune::GeometryTypes::triangle)) + + dofsPerQuad * size_type(gridView_.size(Dune::GeometryTypes::quadrilateral)); + } + case 3: + { + return dofsPerVertex * size_type(gridView_.size(dim)) + + dofsPerEdge * size_type(gridView_.size(dim-1)) + + dofsPerTriangle * size_type(gridView_.size(Dune::GeometryTypes::triangle)) + + dofsPerQuad * size_type(gridView_.size(Dune::GeometryTypes::quadrilateral)) + + dofsPerTetrahedron * size_type(gridView_.size(Dune::GeometryTypes::tetrahedron)) + + dofsPerPyramid * size_type(gridView_.size(Dune::GeometryTypes::pyramid)) + + dofsPerPrism * size_type(gridView_.size(Dune::GeometryTypes::prism)) + + dofsPerHexahedron * size_type(gridView_.size(Dune::GeometryTypes::hexahedron)); + } + } + DUNE_THROW(Dune::NotImplemented, "No size method for " << dim << "d grids available yet!"); + } + + /// Return number of possible values for next position in multi index + size_type size(const SizePrefix prefix) const + { + assert(prefix.size() == 0 || prefix.size() == 1); + return (prefix.size() == 0) ? size() : 0; + } + + /// Get the total dimension of the space spanned by this basis + size_type dimension() const + { + return size(); + } + + /// Get the maximal number of DOFs associated to node for any element + size_type maxNodeSize() const + { + return StaticPower<(k+1),GV::dimension>::power; + } + +protected: + GridView gridView_; + + size_type vertexOffset_; + size_type edgeOffset_; + size_type triangleOffset_; + size_type quadrilateralOffset_; + size_type tetrahedronOffset_; + size_type pyramidOffset_; + size_type prismOffset_; + size_type hexahedronOffset_; +}; + + +template +class GenericLagrangeNode + : public LeafBasisNode +{ + static const int dim = GV::dimension; + static const int maxSize = StaticPower<(k+1),GV::dimension>::power; + + using FiniteElementCache = typename Dune::PQkLocalFiniteElementCache; + +public: + using size_type = std::size_t; + using Element = typename GV::template Codim<0>::Entity; + using FiniteElement = typename FiniteElementCache::FiniteElementType; + + GenericLagrangeNode() = default; + + /// Return current element, throw if unbound + Element const& element() const + { + return *element_; + } + + /// \brief Return the LocalFiniteElement for the element we are bound to + /** + * The LocalFiniteElement implements the corresponding interfaces of the dune-localfunctions module + */ + FiniteElement const& finiteElement() const + { + return *finiteElement_; + } + + /// Bind to element. + void bind(Element const& e) + { + element_ = &e; + finiteElement_ = &(cache_.get(element_->type())); + this->setSize(finiteElement_->size()); + } + +protected: + FiniteElementCache cache_; + FiniteElement const* finiteElement_ = nullptr; + Element const* element_ = nullptr; +}; + + + +template +class GenericLagrangeNodeIndexSet +{ + static const int dim = GV::dimension; + +public: + using size_type = std::size_t; + + /// Type used for global numbering of the basis vectors + using MultiIndex = MI; + + using PreBasis = GenericLagrangePreBasis; + using Node = GenericLagrangeNode; + + GenericLagrangeNodeIndexSet(PreBasis const& preBasis) + : preBasis_(&preBasis) + {} + + /// \brief Bind the view to a grid element + /** + * Having to bind the view to an element before being able to actually access any of its data members + * offers to centralize some expensive setup code in the 'bind' method, which can save a lot of run-time. + */ + void bind(Node const& node) + { + node_ = &node; + } + + /// Unbind the view + void unbind() + { + node_ = nullptr; + } + + /// Size of subtree rooted in this node (element-local) + size_type size() const + { + assert(node_ != nullptr); + return node_->finiteElement().size(); + } + + /// Maps from subtree index set [0..size-1] to a globally unique multi index in global basis + template + Iter indices(Iter it) const + { + assert(node_ != nullptr); + for (size_type i = 0, end = node_->finiteElement().size() ; i < end ; ++it, ++i) + { + Dune::LocalKey localKey = node_->finiteElement().localCoefficients().localKey(i); + const auto& gridIndexSet = preBasis_->gridView().indexSet(); + const auto& element = node_->element(); + + // The dimension of the entity that the current dof is related to + auto dofDim = dim - localKey.codim(); + + // Test for a vertex dof + // The test for k==1 is redundant, but having it here allows the compiler to conclude + // at compile-time that the dofDim==0 case is the only one that will ever happen. + // This leads to measurable speed-up: see + // https://gitlab.dune-project.org/staging/dune-functions/issues/30 + if (k==1 || dofDim==0) { + *it = {{ size_type(gridIndexSet.subIndex(element,localKey.subEntity(),dim)) }}; + continue; + } + + if (dofDim==1) + { // edge dof + if (dim==1) // element dof -- any local numbering is fine + { + *it = {{ preBasis_->edgeOffset_ + + preBasis_->dofsPerEdge * size_type(gridIndexSet.subIndex(element,0,0)) + + localKey.index() }}; + continue; + } + else + { + const auto refElement + = Dune::referenceElement(element.type()); + + // we have to reverse the numbering if the local triangle edge is + // not aligned with the global edge + size_type v0 = gridIndexSet.subIndex(element,refElement.subEntity(localKey.subEntity(),localKey.codim(),0,dim),dim); + size_type v1 = gridIndexSet.subIndex(element,refElement.subEntity(localKey.subEntity(),localKey.codim(),1,dim),dim); + bool flip = (v0 > v1); + *it = {{ (flip) + ? preBasis_->edgeOffset_ + + preBasis_->dofsPerEdge * size_type(gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim())) + + (preBasis_->dofsPerEdge-1)-localKey.index() + : preBasis_->edgeOffset_ + + preBasis_->dofsPerEdge * size_type(gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim())) + + localKey.index() }}; + continue; + } + } + + if (dofDim==2) + { + if (dim==2) // element dof -- any local numbering is fine + { + if (element.type().isTriangle()) + { + const int intLagrangeNodesPerTriangle = (k-1)*(k-2)/2; + *it = {{ preBasis_->triangleOffset_ + + intLagrangeNodesPerTriangle * size_type(gridIndexSet.subIndex(element,0,0)) + + localKey.index() }}; + continue; + } + else if (element.type().isQuadrilateral()) + { + const int intLagrangeNodesPerQuadrilateral = (k-1)*(k-1); + *it = {{ preBasis_->quadrilateralOffset_ + + intLagrangeNodesPerQuadrilateral * size_type(gridIndexSet.subIndex(element,0,0)) + + localKey.index() }}; + continue; + } + else + DUNE_THROW(Dune::NotImplemented, "2d elements have to be triangles or quadrilaterals"); + } + else + { + const auto refElement + = Dune::referenceElement(element.type()); + + if (k>3) + DUNE_THROW(Dune::NotImplemented, "GenericLagrangeNodalBasis for 3D grids is only implemented if k<=3"); + + if (k==3 and !refElement.type(localKey.subEntity(), localKey.codim()).isTriangle()) + DUNE_THROW(Dune::NotImplemented, "GenericLagrangeNodalBasis for 3D grids with k==3 is only implemented if the grid is a simplex grid"); + + *it = {{ preBasis_->triangleOffset_ + + size_type(gridIndexSet.subIndex(element,localKey.subEntity(),localKey.codim())) }}; + continue; + } + } + + if (dofDim==3) + { + if (dim==3) // element dof -- any local numbering is fine + { + if (element.type().isTetrahedron()) + { + *it = {{ preBasis_->tetrahedronOffset_ + + PreBasis::dofsPerTetrahedron*((size_type)gridIndexSet.subIndex(element,0,0)) + + localKey.index() }}; + continue; + } + else if (element.type().isHexahedron()) + { + *it = {{ preBasis_->hexahedronOffset_ + + PreBasis::dofsPerHexahedron*((size_type)gridIndexSet.subIndex(element,0,0)) + + localKey.index() }}; + continue; + } + else if (element.type().isPrism()) + { + *it = {{ preBasis_->prismOffset_ + + PreBasis::dofsPerPrism*((size_type)gridIndexSet.subIndex(element,0,0)) + + localKey.index() }}; + continue; + } + else if (element.type().isPyramid()) + { + *it = {{ preBasis_->pyramidOffset_ + + PreBasis::dofsPerPyramid*((size_type)gridIndexSet.subIndex(element,0,0)) + + localKey.index() }}; + continue; + } + else + DUNE_THROW(Dune::NotImplemented, "3d elements have to be tetrahedra, hexahedra, prisms, or pyramids"); + } else + DUNE_THROW(Dune::NotImplemented, "Grids of dimension larger than 3 are no supported"); + } + DUNE_THROW(Dune::NotImplemented, "Grid contains elements not supported for the GenericLagrangeNodalBasis"); + } + return it; + } + +protected: + PreBasis const* preBasis_ = nullptr; + Node const* node_ = nullptr; +}; + + +namespace BasisFactory { +namespace Imp { + +template +class GenericLagrangePreBasisFactory +{ +public: + static const std::size_t requiredMultiIndexSize = 1; + + template + auto makePreBasis(GridView const& gridView) const + { + return GenericLagrangePreBasis(gridView); + } +}; + +} // end namespace BasisFactory::Imp + + +/// \brief Create a pre-basis factory that can create a GenericLagrange pre-basis +/** + * \ingroup FunctionSpaceBasesImplementations + * + * \tparam k The polynomial order of ansatz functions + * \tparam R The range type of the local basis + */ +template +auto lagrange() +{ + return Imp::GenericLagrangePreBasisFactory(); +} + +} // end namespace BasisFactory + + +/// \brief Nodal basis of a scalar k-th-order Lagrangean finite element space with generic range type +/** + * \ingroup FunctionSpaceBasesImplementations + * + * \note This only works for certain grids. The following restrictions hold + * - If k is no larger than 2, then the grids can have any dimension + * - If k is larger than 3 then the grid must be two-dimensional + * - If k is 3, then the grid can be 3d *if* it is a simplex grid + * + * All arguments passed to the constructor will be forwarded to the constructor + * of GenericLagrangePreBasis. + * + * \tparam GV The GridView that the space is defined on + * \tparam k The order of the basis + * \tparam R The range type of the local basis + */ +template +using GenericLagrangeBasis = DefaultGlobalBasis> >; + +} // end namespace Functions +} // end namespace Dune -- GitLab From c485820cf147c4a0e035266f3978b499c826055c Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Tue, 3 Sep 2019 13:12:40 +0200 Subject: [PATCH 2/3] added generic lagrange basis to ProblemStatTest --- test/ProblemStatTest.cpp | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/test/ProblemStatTest.cpp b/test/ProblemStatTest.cpp index 8c6b36db..450fc1fc 100644 --- a/test/ProblemStatTest.cpp +++ b/test/ProblemStatTest.cpp @@ -2,20 +2,25 @@ #include #include #include +#include using namespace AMDiS; -using Grid = Dune::YaspGrid<2>; - -template -struct Param - : DefaultBasisCreator, T> -{}; - template void test() { - ProblemStat> prob("ellipt"); + // use T as coordinate type + using Grid = Dune::YaspGrid<2, Dune::EquidistantCoordinates>; + + // use T as range type for basis + using namespace Dune::Function::BasisFactory; + auto basis = makeBasis(grid.leafGridView(), power<1>(lagrange<1,T>(), flatLexicographic())); + using Basis = decltype(basis); + + // use T as coefficient type + using Param = DefaultProblemTraits; + + ProblemStat prob("ellipt"); prob.initialize(INIT_ALL); prob.boundaryManager()->setBoxBoundary({1,1,2,2}); -- GitLab From 4192b5f90242f18042674cd4f8a265eb5ff1a993 Mon Sep 17 00:00:00 2001 From: Simon Praetorius Date: Tue, 3 Sep 2019 13:35:58 +0200 Subject: [PATCH 3/3] added dune 2.6 implementation for GenericLagrangeBasis --- src/amdis/functions/GenericLagrangeBasis.hpp | 74 ++++++++++++++------ test/ProblemStatTest.cpp | 7 +- 2 files changed, 58 insertions(+), 23 deletions(-) diff --git a/src/amdis/functions/GenericLagrangeBasis.hpp b/src/amdis/functions/GenericLagrangeBasis.hpp index 9601a84e..986060af 100644 --- a/src/amdis/functions/GenericLagrangeBasis.hpp +++ b/src/amdis/functions/GenericLagrangeBasis.hpp @@ -9,10 +9,6 @@ #include #include -#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7) -#error "This implementation of GenericLagrange basis requires dune 2.7" -#endif - namespace Dune { namespace Functions { @@ -33,17 +29,24 @@ namespace Functions { **/ // forward declarations + +#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7) +template +class GenericLagrangeNode; + +template +class GenericLagrangeNodeIndexSet; +#else template class GenericLagrangeNode; template class GenericLagrangeNodeIndexSet; +#endif template class GenericLagrangePreBasis; - - /** * \brief A pre-basis for a PQ-lagrange bases with given order and range type * @@ -72,8 +75,13 @@ public: using size_type = std::size_t; private: +#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7) + template + friend class GenericLagrangeNodeIndexSet; +#else template friend class GenericLagrangeNodeIndexSet; +#endif // Precompute the number of dofs per entity type const static size_type dofsPerVertex = @@ -94,11 +102,18 @@ private: k == 0 ? (dim == 3 ? 1 : 0) : (k-2)*(k-1)*(2*k-3)/6; public: + +#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7) + template + using Node = GenericLagrangeNode; + template + using IndexSet = LagrangeNodeIndexSet; +#else /// Template mapping root tree path to type of created tree node using Node = GenericLagrangeNode; - /// Template mapping root tree path to type of created tree node index set using IndexSet = GenericLagrangeNodeIndexSet; +#endif /// Type used for global numbering of the basis vectors using MultiIndex = MI; @@ -142,22 +157,19 @@ public: gridView_ = gridView; } +#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7) + template + Node node(const TP& tp) const { return Node{tp}; } + + template + IndexSet indexSet() const { return IndexSet{*this}; } +#else /// Create tree node - Node makeNode() const - { - return Node{}; - } + Node makeNode() const { return Node{}; } - /** - * \brief Create tree node index set - * - * Create an index set suitable for the tree node obtained - * by makeNode(). - */ - IndexSet makeIndexSet() const - { - return IndexSet{*this}; - } + /// Create tree node index set + IndexSet makeIndexSet() const { return IndexSet{*this}; } +#endif /// Same as size(prefix) with empty prefix size_type size() const @@ -223,9 +235,15 @@ protected: }; +#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7) +template +class GenericLagrangeNode + : public LeafBasisNode +#else template class GenericLagrangeNode : public LeafBasisNode +#endif { static const int dim = GV::dimension; static const int maxSize = StaticPower<(k+1),GV::dimension>::power; @@ -237,7 +255,12 @@ public: using Element = typename GV::template Codim<0>::Entity; using FiniteElement = typename FiniteElementCache::FiniteElementType; +#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7) + using TreePath = TP; + GenericLagrangeNode(const TreePath& treePath) : LeafBasisNode(treePath) {} +#else GenericLagrangeNode() = default; +#endif /// Return current element, throw if unbound Element const& element() const @@ -270,7 +293,11 @@ protected: +#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7) +template +#else template +#endif class GenericLagrangeNodeIndexSet { static const int dim = GV::dimension; @@ -282,7 +309,12 @@ public: using MultiIndex = MI; using PreBasis = GenericLagrangePreBasis; + +#if DUNE_VERSION_LT(DUNE_FUNCTIONS,2,7) + using Node = typename PreBasis::template Node; +#else using Node = GenericLagrangeNode; +#endif GenericLagrangeNodeIndexSet(PreBasis const& preBasis) : preBasis_(&preBasis) diff --git a/test/ProblemStatTest.cpp b/test/ProblemStatTest.cpp index 450fc1fc..f2ecc682 100644 --- a/test/ProblemStatTest.cpp +++ b/test/ProblemStatTest.cpp @@ -4,6 +4,8 @@ #include #include +#include + using namespace AMDiS; template @@ -11,16 +13,17 @@ void test() { // use T as coordinate type using Grid = Dune::YaspGrid<2, Dune::EquidistantCoordinates>; + Grid grid({T(1), T(1)}, {8,8}); // use T as range type for basis - using namespace Dune::Function::BasisFactory; + using namespace Dune::Functions::BasisFactory; auto basis = makeBasis(grid.leafGridView(), power<1>(lagrange<1,T>(), flatLexicographic())); using Basis = decltype(basis); // use T as coefficient type using Param = DefaultProblemTraits; - ProblemStat prob("ellipt"); + ProblemStat prob("ellipt", grid, basis); prob.initialize(INIT_ALL); prob.boundaryManager()->setBoxBoundary({1,1,2,2}); -- GitLab