Loading...
Searching...
No Matches
morleybasis.hh
Go to the documentation of this file.
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_MORLEYBASIS_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_MORLEYBASIS_HH
#include <algorithm>
#include <type_traits>
#include <vector>
#include <array>
#include <bitset>
#include <dune/common/boundschecking.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/common/rangegenerators.hh>
#include <dune/localfunctions/common/localbasis.hh>
#include <dune/localfunctions/common/localfiniteelementtraits.hh>
#include <dune/localfunctions/common/localkey.hh>
#include <dune/functions/analyticfunctions/monomialset.hh>
#include <dune/functions/common/densevectorview.hh>
#include <dune/functions/common/squeezetensor.hh>
#include <dune/functions/functionspacebases/cubichermitebasis.hh>
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/functionaldescriptor.hh>
#include <dune/functions/functionspacebases/leafprebasismappermixin.hh>
#include <dune/functions/functionspacebases/nodes.hh>
#include <dune/functions/functionspacebases/transformedfiniteelementmixin.hh>
namespace Dune::Functions
{
template<class GV, class R>
class MorleyPreBasis;
template<class GV, class R = double>
using MorleyBasis = DefaultGlobalBasis<MorleyPreBasis<GV, R> >;
namespace Impl
{
template<class D, class R>
class MorleyReferenceLocalBasis
{
public:
static constexpr int dim = 2;
using Traits = H2LocalBasisTraits<D, dim, FieldVector<D, dim>, R, 1, FieldVector<R, 1>,
FieldMatrix<R, 1, dim>, FieldMatrix<R, dim, dim>>;
private:
static constexpr auto getMorleyCoefficients()
{
// Define std::sqrt(2.) manually in double precision,
// because std::sqrt is not constexpr before C++26.
D sqrt2 = 0.5 * 1.414213562373095;
return Dune::FieldMatrix<D, 6,6>({{1, -1, -1, 0, 2, 0},
{0, 0.5, 0.5, 0.5, -1, -0.5},
{0, 0.5, 0.5, -0.5, -1, 0.5},
{0, 0, 1, 0, 0, -1},
{0, -1., 0, 1., 0, 0},
{0, sqrt2, sqrt2, -sqrt2, -2. * sqrt2, -sqrt2}});
}
static constexpr auto referenceBasisCoefficients = getMorleyCoefficients();
MonomialSet<typename Traits::RangeFieldType, dim, 2> monomials;
public:
static constexpr unsigned int size()
{
return 6;
}
static constexpr unsigned int order()
{
return 2;
}
void evaluateFunction(const typename Traits::DomainType &in,
std::vector<typename Traits::RangeType> &out) const
{
out.resize(size());
auto monomialValues = monomials(in);
multiplyWithCoefficentMatrix(referenceBasisCoefficients, monomialValues, out);
}
void evaluateJacobian(const typename Traits::DomainType &in,
std::vector<typename Traits::JacobianType> &out) const
{
out.resize(size());
auto monomialValues = derivative(monomials)(in);
multiplyWithCoefficentMatrix(referenceBasisCoefficients, monomialValues, out);
}
void evaluateHessian(const typename Traits::DomainType &in,
std::vector<typename Traits::HessianType> &out) const
{
out.resize(size());
auto monomialValues = derivative(derivative(monomials))(in);
multiplyWithCoefficentMatrix(referenceBasisCoefficients, monomialValues, out);
}
void partial(std::array<unsigned int, dim> order, const typename Traits::DomainType &in,
std::vector<typename Traits::RangeType> &out) const
{
out.resize(size());
auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
if (totalOrder == 0)
evaluateFunction(in, out);
else if (totalOrder == 1)
{
evaluateJacobian(in,jacobiansBuffer_);
std::size_t which = std::max_element(order.begin(), order.end()) - order.begin();
for (auto i : Dune::range(size()))
out[i] = jacobiansBuffer_[i][0][which];
}
else if (totalOrder == 2)
{
evaluateHessian(in, hessianBuffer_);
std::size_t first, second;
first = std::max_element(order.begin(), order.end()) - order.begin();
if (order[first] == 2)
second = first;
else
{
order[first] = 0;
second = std::max_element(order.begin(), order.end()) - order.begin();
}
for (auto i : Dune::range(size()))
out[i] = hessianBuffer_[i][first][second];
}
else
DUNE_THROW(RangeError, "partial() not implemented for given order");
}
private:
mutable std::vector<typename Traits::JacobianType> jacobiansBuffer_;
mutable std::vector<typename Traits::HessianType> hessianBuffer_;
};
class MorleyLocalCoefficients
{
public:
using size_type = std::size_t;
MorleyLocalCoefficients()
{
for (size_type i = 0; i < 3; ++i)
{
localKeys_[i] = LocalKey(i, 2, 0); // vertex dofs
localKeys_[3 + i] = LocalKey(i, 1, 0); // edge dofs
}
}
static constexpr size_type size()
{
return 6;
}
constexpr LocalKey const &localKey(size_type i) const
{
return localKeys_[i];
}
private:
std::array<LocalKey, 6> localKeys_;
};
template<class D>
class MorleyLocalInterpolation
{
using size_type = std::size_t;
using FunctionalDescriptor = Dune::Functions::Impl::FunctionalDescriptor<2>;
public:
MorleyLocalInterpolation()
{
descriptors_[0] = FunctionalDescriptor();
descriptors_[1] = FunctionalDescriptor();
descriptors_[2] = FunctionalDescriptor();
descriptors_[3] = FunctionalDescriptor(1);
descriptors_[4] = FunctionalDescriptor(1);
descriptors_[5] = FunctionalDescriptor(1);
}
template <class Element>
void bind(const Element &element, const std::bitset<3> &edgeOrientation)
{
auto geometry = element.geometry();
// get global Normals and midpoints
auto refElement = Dune::referenceElement<double, 2>(geometry.type());
for (std::size_t i = 0; i < 3; ++i)
{
localVertices_[i] = refElement.position(i, 2);
localMidpoints_[i] = refElement.position(i, 1);
std::size_t lower = (i == 2) ? 1 : 0;
std::size_t upper = (i == 0) ? 1 : 2;
auto edge = geometry.global(refElement.position(upper, 2))
- geometry.global(refElement.position(lower, 2));
// normalize and orient
edge /= edge.two_norm() * (edgeOrientation[i] ? -1. : 1.);
// Rotation by pi/2. Note that Kirby rotates by 3*pi/2
globalNormals_[i] = {-edge[1], edge[0]};
}
}
template <class F, class C>
void interpolate(const F &f, std::vector<C> &out) const
{
out.resize(6);
auto&& df = derivative(f);
for (size_type i = 0; i < 3; ++i)
{
out[i] = f(localVertices_[i]);
out[3 + i] = squeezeTensor(df(localMidpoints_[i])).dot(globalNormals_[i]);
}
}
const FunctionalDescriptor& functionalDescriptor(size_type i) const
{
return descriptors_[i];
}
protected:
std::array<Dune::FieldVector<D, 2>, 3> globalNormals_;
std::array<Dune::FieldVector<D, 2>, 3> localMidpoints_;
std::array<Dune::FieldVector<D, 2>, 3> localVertices_;
std::array<FunctionalDescriptor, 6> descriptors_;
};
template<class D, class R>
using MorleyLocalBasisTraits = typename Impl::MorleyReferenceLocalBasis<D, R>::Traits;
template<class D, class R>
class MorleyLocalFiniteElement
: public Impl::TransformedFiniteElementMixin<MorleyLocalFiniteElement<D,R>, MorleyLocalBasisTraits<D, R>>
{
using Base = Impl::TransformedFiniteElementMixin< MorleyLocalFiniteElement<D,R>, MorleyLocalBasisTraits<D, R>>;
friend class Impl::TransformedLocalBasis<MorleyLocalFiniteElement<D,R>, MorleyLocalBasisTraits<D, R>>;
static constexpr int dim = 2;
public:
using size_type = std::size_t;
using Traits = LocalFiniteElementTraits<
Impl::TransformedLocalBasis<MorleyLocalFiniteElement<D,R>, MorleyLocalBasisTraits<D, R>>,
Impl::MorleyLocalCoefficients,
Impl::MorleyLocalInterpolation<D>>;
const typename Traits::LocalCoefficientsType &localCoefficients() const
{
return coefficients_;
}
const typename Traits::LocalInterpolationType &localInterpolation() const
{
return interpolation_;
}
static constexpr GeometryType type()
{
return GeometryTypes::simplex(dim);
}
static constexpr size_type size()
{
return 6;
}
template<class Element>
void bind(std::bitset<3> const& data, Element const &e)
{
edgeOrientation_ = data;
fillMatrix(e.geometry());
interpolation_.bind(e, edgeOrientation_);
}
protected:
Impl::MorleyReferenceLocalBasis<D, R> const& referenceLocalBasis() const
{
return basis_;
}
template<class InputValues, class OutputValues>
void transform(InputValues const& inValues, OutputValues& outValues) const
{
// Here we cannot directly use
// mat_.mv(inValues, outValues);
// because mv expects the DenseVector interface.
auto inValuesDenseVector = Impl::DenseVectorView(inValues);
auto outValuesDenseVector = Impl::DenseVectorView(outValues);
mat_.mv(inValuesDenseVector, outValuesDenseVector);
}
private:
template<class Geometry>
void fillMatrix(Geometry const &geometry)
{
std::array<R, 3> B_11;
std::array<R, 3> B_12;
std::array<R, 3> l_inv;
std::array<Dune::FieldVector<R, 2>, 3> referenceTangents;
std::array<Dune::FieldVector<R, 2>, 3> globalTangents;
// By default, edges point from the vertex with the smaller index
// to the vertex with the larger index.
// get local and global Tangents
auto refElement = Dune::referenceElement<double, 2>(geometry.type());
auto x = refElement.position(0,0);
for (std::size_t i = 0; i < 3; ++i)
{
std::size_t lower = (i == 2) ? 1 : 0;
std::size_t upper = (i == 0) ? 1 : 2;
auto edge = refElement.position(upper, 2) - refElement.position(lower, 2);
referenceTangents[i] = edge / edge.two_norm();
auto globalEdge = geometry.global(refElement.position(upper, 2))
- geometry.global(refElement.position(lower, 2));
l_inv[i] = 1. / globalEdge.two_norm();
globalTangents[i] = globalEdge * l_inv[i];
}
auto jacobianTransposed = geometry.jacobianTransposed(x);
for (std::size_t i = 0; i < 3; ++i)
{
B_11[i] = -referenceTangents[i][1]
*(-globalTangents[i][1] * jacobianTransposed[0][0]
+ globalTangents[i][0] * jacobianTransposed[0][1])
+ referenceTangents[i][0]
*(-globalTangents[i][1] * jacobianTransposed[1][0]
+ globalTangents[i][0] * jacobianTransposed[1][1]);
B_12[i] = -referenceTangents[i][1]
*(globalTangents[i][0] * jacobianTransposed[0][0]
+ globalTangents[i][1] * jacobianTransposed[0][1])
+ referenceTangents[i][0]
*(globalTangents[i][0] * jacobianTransposed[1][0]
+ globalTangents[i][1] * jacobianTransposed[1][1]);
}
// Actually setup matrix
int sign = -1;
for (std::size_t i = 0; i < 3; ++i)
{
mat_[i][i] = 1.;
for (std::size_t j = 0; j < 3; ++j)
{
if (j != (2 - i)) // dune specific edge order
{
mat_[j][3 + i] = sign * B_12[i] * l_inv[i];
sign *= -1;
}
}
mat_[3 + i][3 + i] = (edgeOrientation_[i] ? -1. : 1.) * B_11[i];
}
}
// a finite element consists of a basis, coeffiecents and an interpolation
typename Impl::MorleyReferenceLocalBasis<D, R> basis_;
typename Traits::LocalCoefficientsType coefficients_;
typename Traits::LocalInterpolationType interpolation_;
// This is the matrix M in Kirbys paper
Dune::FieldMatrix<R, 6, 6>mat_;
// the local state, i.e. a collection of global information restricted to this element
std::bitset<3> edgeOrientation_;
};
} // end namespace Impl
// *****************************************************************************
// This is the reusable part of the basis. It contains
//
// MorleyPreBasis
// MorleyNode
//
// The pre-basis allows to create the others and is the owner of possible shared
// state. These components do _not_ depend on the global basis and local view
// and can be used without a global basis.
// *****************************************************************************
template<class GV, class R>
class MorleyNode
: public LeafBasisNode
{
using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GV>;
public:
using size_type = std::size_t;
using Element = typename GV::template Codim<0>::Entity;
using FiniteElement = typename Impl::MorleyLocalFiniteElement<typename GV::ctype, R>;
MorleyNode(Mapper const& m, std::vector<std::bitset<3>> const& data)
: mapper_(&m), data_(&data)
{}
Element const &element() const
{
return *element_;
}
FiniteElement const &finiteElement() const
{
return finiteElement_;
}
void bind(Element const &e)
{
element_ = &e;
finiteElement_.bind((*data_)[mapper_->index(e)], *element_);
this->setSize(finiteElement_.size());
}
unsigned int order() const
{
return finiteElement_.localBasis().order();
}
protected:
FiniteElement finiteElement_;
Element const* element_;
Mapper const* mapper_;
std::vector<std::bitset<3>> const* data_;
};
template<class GV, class R>
class MorleyPreBasis
: public LeafPreBasisMapperMixin<GV>
{
using Base = LeafPreBasisMapperMixin<GV>;
using SubEntityMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GV>;
using Element = typename GV::template Codim<0>::Entity;
using D = typename GV::ctype;
static const std::size_t dim = GV::dimension;
// helper methods to assign each subentity the number of dofs. Used by the LeafPreBasisMapperMixin.
static constexpr auto morleyMapperLayout(Dune::GeometryType type, int gridDim)
{
assert(gridDim == 2);
if (type.isVertex())
return 1; // one evaluation dof per vertex
if (type.isLine())
return 1;
if ((type.isTriangle()) )
return 0;
else
return 0;
}
public:
using GridView = GV;
using size_type = std::size_t;
using Node = MorleyNode<GridView, R>;
MorleyPreBasis(const GV &gv)
: Base(gv, morleyMapperLayout)
, mapper_({gv, mcmgElementLayout()})
{
data_ = Impl::computeEdgeOrientations(mapper_);
}
void update(GridView const &gv)
{
Base::update(gv);
data_ = Impl::computeEdgeOrientations(mapper_);
}
Node makeNode() const
{
return Node{mapper_, data_};
}
protected:
SubEntityMapper mapper_;
std::vector<std::bitset<3>> data_;
}; // class MorleyPreBasis
namespace BasisFactory
{
template<class R = double>
auto morley()
{
return [=](auto const &gridView) {
return MorleyPreBasis<std::decay_t<decltype(gridView)>, R>(gridView);
};
}
} // end namespace BasisFactory
} // end namespace Dune::Functions
#endif
nodes.hh
leafprebasismappermixin.hh
functionaldescriptor.hh
cubichermitebasis.hh
This file provides an implementation of the cubic Hermite finite element in 1 to 3 dimensions.
defaultglobalbasis.hh
transformedfiniteelementmixin.hh
monomialset.hh
densevectorview.hh
squeezetensor.hh
Dune::Functions::derivative
TrigonometricFunction< K, -cosFactor, sinFactor > derivative(const TrigonometricFunction< K, sinFactor, cosFactor > &f)
Obtain derivative of TrigonometricFunction function.
Definition trigonometricfunction.hh:43
Dune::Functions::MorleyBasis
DefaultGlobalBasis< MorleyPreBasis< GV, R > > MorleyBasis
Nodal basis of a scalar quadratic Morley finite element space.
Definition morleybasis.hh:76
Dune::Functions
Definition monomialset.hh:19
Dune::Functions::interpolate
void interpolate(const B &basis, C &&coeff, const F &f, const BV &bv, const NTRE &nodeToRangeEntry)
Interpolate given function in discrete function space.
Definition interpolate.hh:205
Dune::Functions::BasisFactory
Definition argyrisbasis.hh:926
Dune::Functions::BasisFactory::morley
auto morley()
construct a PreBasisFactory for the Morley Finite Element
Definition morleybasis.hh:672
Dune::Functions::MonomialSet
Function, which evaluates all monomials up to degree maxDegree in a given coordinate.
Definition monomialset.hh:64
Dune::Functions::H2LocalBasisTraits
Definition cubichermitebasis.hh:84
Dune::Functions::DefaultGlobalBasis
Global basis for given pre-basis.
Definition defaultglobalbasis.hh:53
Dune::Functions::LeafPreBasisMapperMixin::update
void update(const GridView &gv)
Update the stored GridView.
Definition leafprebasismappermixin.hh:101
Dune::Functions::LeafPreBasisMapperMixin::LeafPreBasisMapperMixin
LeafPreBasisMapperMixin(const GridView &gv, Dune::MCMGLayout layout)
Construct from GridView and local DOF layout.
Definition leafprebasismappermixin.hh:74
Dune::Functions::MorleyPreBasis
A pre-basis for a Morleybasis.
Definition morleybasis.hh:601
Dune::Functions::MorleyPreBasis::update
void update(GridView const &gv)
Update the stored grid view, to be called if the grid has changed.
Definition morleybasis.hh:641
Dune::Functions::MorleyPreBasis::mapper_
SubEntityMapper mapper_
Definition morleybasis.hh:657
Dune::Functions::MorleyPreBasis::makeNode
Node makeNode() const
Create tree node.
Definition morleybasis.hh:650
Dune::Functions::MorleyPreBasis::data_
std::vector< std::bitset< 3 > > data_
Definition morleybasis.hh:658
Dune::Functions::MorleyPreBasis::size_type
std::size_t size_type
Type used for indices and size information.
Definition morleybasis.hh:627
Dune::Functions::MorleyPreBasis::GridView
GV GridView
The grid view that the FE basis is defined on.
Definition morleybasis.hh:624
Dune::Functions::MorleyPreBasis::Node
MorleyNode< GridView, R > Node
Template mapping root tree path to type of created tree node.
Definition morleybasis.hh:630
Dune::Functions::MorleyPreBasis::MorleyPreBasis
MorleyPreBasis(const GV &gv)
Constructor for a given grid view object.
Definition morleybasis.hh:633
Dune::Functions::MorleyNode
Definition morleybasis.hh:539
Dune::Functions::MorleyNode::element
Element const & element() const
Return current element, throw if unbound.
Definition morleybasis.hh:552
Dune::Functions::MorleyNode::size_type
std::size_t size_type
Definition morleybasis.hh:542
Dune::Functions::MorleyNode::finiteElement
FiniteElement const & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition morleybasis.hh:562
Dune::Functions::MorleyNode< GridView, R >::element_
Element const * element_
Definition morleybasis.hh:583
Dune::Functions::MorleyNode::order
unsigned int order() const
The order of the local basis.
Definition morleybasis.hh:576
Dune::Functions::MorleyNode< GridView, R >::data_
std::vector< std::bitset< 3 > > const * data_
Definition morleybasis.hh:585
Dune::Functions::MorleyNode< GridView, R >::mapper_
Mapper const * mapper_
Definition morleybasis.hh:584
Dune::Functions::MorleyNode::Element
typename GV::template Codim< 0 >::Entity Element
Definition morleybasis.hh:543
Dune::Functions::MorleyNode::bind
void bind(Element const &e)
Bind to element.
Definition morleybasis.hh:568
Dune::Functions::MorleyNode::MorleyNode
MorleyNode(Mapper const &m, std::vector< std::bitset< 3 > > const &data)
Definition morleybasis.hh:547
Dune::Functions::MorleyNode::FiniteElement
typename Impl::MorleyLocalFiniteElement< typename GV::ctype, R > FiniteElement
Definition morleybasis.hh:544
Dune::Functions::MorleyNode< GridView, R >::finiteElement_
FiniteElement finiteElement_
Definition morleybasis.hh:582
Dune::Functions::BasisNodeMixin::setSize
void setSize(const size_type size)
Definition nodes.hh:197
Dune::Functions::LeafBasisNode
Definition nodes.hh:218