Loading...
Searching...
No Matches
monomialset.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_ANALYTICFUNCTIONS_MONOMIALSET_HH
#define DUNE_FUNCTIONS_ANALYTICFUNCTIONS_MONOMIALSET_HH
#include <array>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/math.hh>
#include <dune/common/rangeutilities.hh>
namespace Dune::Functions {
namespace Impl {
// Computes `y = [1, x, x^2, x^3, ..., x^maxDegree]` with `x` a single coordinate
template<int maxDegree, class DomainFieldType, class RangeType>
void computePowers(const DomainFieldType& x, RangeType& y)
{
if constexpr(maxDegree >= 0)
{
y[0] = 1;
for(auto k : Dune::range(maxDegree))
y[k+1] = y[k]*x;
}
}
} // namespace Impl
template<class RangeFieldType, int dimension, int maxDegree>
struct MonomialSet
{
static constexpr int dim = dimension;
static constexpr int size = Dune::binomial(maxDegree + dim, dim);
static_assert(1 <= dim && dim <= 3, "MonomialSet is specialized for dimensions 1,2,3 only.");
static constexpr std::array<std::array<std::size_t,dim>,size> exponents();
template<class DomainFieldType>
constexpr Dune::FieldVector<RangeFieldType,size> operator()(const Dune::FieldVector<DomainFieldType,dim>& x) const;
struct Derivative
{
template<class DomainFieldType>
constexpr Dune::FieldMatrix<RangeFieldType,size, dim> operator()(const Dune::FieldVector<DomainFieldType,dim>& x) const;
struct Hessian
{
template<class DomainFieldType>
constexpr std::array<FieldMatrix<RangeFieldType,dim, dim>, size> operator()(const Dune::FieldVector<DomainFieldType,dim>& x) const;
};
constexpr friend auto derivative(const Derivative & d)
{
return Hessian{};
}
};
constexpr friend auto derivative(const MonomialSet& m)
{
return Derivative{};
}
};
#ifndef DOXYGEN
// Specialization for dim = 1
template<class RangeFieldType, int maxDegree>
struct MonomialSet<RangeFieldType, 1, maxDegree>
{
static constexpr int dim = 1;
static constexpr int size = maxDegree+1;
static constexpr auto exponents()
{
auto p = std::array<std::array<std::size_t,1>,size>{};
for(auto k : Dune::range(size))
p[k][0] = k;
return p;
}
template<class DomainFieldType>
constexpr auto operator()(const Dune::FieldVector<DomainFieldType,1>& x) const
{
auto y = Dune::FieldVector<RangeFieldType,size>{};
Impl::computePowers<maxDegree>(x[0], y);
return y;
}
struct Derivative
{
template<class DomainFieldType>
constexpr auto operator()(const Dune::FieldVector<DomainFieldType,1>& x) const
{
auto xPowers = Dune::FieldVector<RangeFieldType,size>{};
Impl::computePowers<maxDegree-1>(x[0], xPowers);
auto y = Dune::FieldMatrix<RangeFieldType,size,1>{};
for(auto degree : Dune::range(1, maxDegree+1))
y[degree][0] = degree*xPowers[degree-1];
return y;
}
struct Hessian
{
template<class DomainFieldType>
constexpr auto operator()(const Dune::FieldVector<DomainFieldType,1>& x) const
{
auto xPowers = std::array<RangeFieldType,maxDegree+1>{};
Impl::computePowers<maxDegree-2>(x[0],xPowers);
auto y = std::array<Dune::FieldMatrix<RangeFieldType,1,1>,size>{};
for(auto degree : Dune::range(maxDegree+1))
if (degree-1 > 0)
y[degree][0][0] = degree*(degree-1)*xPowers[degree-2];
return y;
}
};
constexpr friend auto derivative(const Derivative& d)
{
return Hessian{};
}
};
constexpr friend auto derivative(const MonomialSet& m)
{
return Derivative{};
}
};
// Specialization for dim = 2
template<class RangeFieldType, int maxDegree>
struct MonomialSet<RangeFieldType, 2, maxDegree>
{
static constexpr int dim = 2;
static constexpr int size = (maxDegree+1)*(maxDegree+2)/2;
static constexpr auto exponents()
{
auto p = std::array<std::array<std::size_t,2>,size>{};
std::size_t index = 0;
for(auto degree : Dune::range(maxDegree+1))
{
for(auto k : Dune::range(degree+1))
{
p[index][0] = degree-k;
p[index][1] = k;
++index;
}
}
return p;
}
template<class DomainFieldType>
constexpr auto operator()(const Dune::FieldVector<DomainFieldType,2>& x) const
{
auto xPowers = std::array<Dune::FieldVector<RangeFieldType,size>,dim>{};
for(auto j : Dune::range(dim))
Impl::computePowers<maxDegree>(x[j], xPowers[j]);
auto y = Dune::FieldVector<RangeFieldType,size>{};
std::size_t index = 0;
for(auto degree : Dune::range(maxDegree+1))
{
for(auto k : Dune::range(degree+1))
{
y[index] = xPowers[0][degree-k]*xPowers[1][k];
++index;
}
}
return y;
}
struct Derivative
{
template<class DomainFieldType>
constexpr auto operator()(const Dune::FieldVector<DomainFieldType,2>& x) const
{
auto xPowers = std::array<Dune::FieldVector<RangeFieldType,size>,dim>{};
for(auto j : Dune::range(dim))
Impl::computePowers<maxDegree-1>(x[j], xPowers[j]);
auto y = Dune::FieldMatrix<RangeFieldType,size,2>{};
std::size_t index = 0;
for(auto degree : Dune::range(maxDegree+1))
{
for(auto k : Dune::range(degree+1))
{
if (degree-k > 0)
y[index][0] = (degree-k)*xPowers[0][degree-k-1]*xPowers[1][k];
if (k > 0)
y[index][1] = k*xPowers[0][degree-k]*xPowers[1][k-1];
++index;
}
}
return y;
}
struct Hessian
{
template<class DomainFieldType>
constexpr auto operator()(const Dune::FieldVector<DomainFieldType,2>& x) const
{
auto xPowers = std::array<std::array<RangeFieldType,maxDegree+1>,dim>{};
for(auto j : Dune::range(dim))
Impl::computePowers<maxDegree-2>(x[j], xPowers[j]);
auto y = std::array<Dune::FieldMatrix<RangeFieldType,2,2>,size>{};
std::size_t index = 0;
for(auto degree : Dune::range(maxDegree+1))
{
for(auto k : Dune::range(degree+1))
{
if (degree-k > 1)
y[index][0][0] = (degree-k-1)*(degree-k)*xPowers[0][degree-k-2]*xPowers[1][k];
if (k > 0 and degree-k > 0){
auto mixed = k*(degree-k)*xPowers[0][degree-k-1]*xPowers[1][k-1];
y[index][0][1]= mixed;
y[index][1][0]= mixed;
}
if (k > 1)
y[index][1][1] = k*(k-1)*xPowers[0][degree-k]*xPowers[1][k-2];
++index;
}
}
return y;
}
};
constexpr friend auto derivative(const Derivative & d)
{
return Hessian{};
}
};
constexpr friend auto derivative(const MonomialSet& m)
{
return Derivative{};
}
};
// Specialization for dim = 3
template<class RangeFieldType, int maxDegree>
struct MonomialSet<RangeFieldType, 3, maxDegree>
{
static constexpr int dim = 3;
static constexpr int size = Dune::binomial(maxDegree + dim, dim);
static constexpr auto exponents()
{
auto p = std::array<std::array<std::size_t,3>,size>{};
std::size_t index = 0;
for(auto degree : Dune::range(maxDegree+1))
{
for(auto k : Dune::range(degree+1))
{
for (auto l : Dune::range(degree-k+1))
{
p[index][0] = degree-k-l;
p[index][1] = l;
p[index][2] = k;
++index;
}
}
}
return p;
}
template<class DomainFieldType>
constexpr auto operator()(const Dune::FieldVector<DomainFieldType,3>& x) const
{
auto xPowers = std::array<std::array<RangeFieldType,maxDegree+1>,dim>{};
for(auto j : Dune::range(dim))
Impl::computePowers<maxDegree>(x[j], xPowers[j]);
auto y = Dune::FieldVector<RangeFieldType,size>{};
std::size_t index = 0;
for(auto degree : Dune::range(maxDegree+1))
{
for(auto k : Dune::range(degree+1))
{
for (auto l : Dune::range(degree-k+1))
{
y[index] = xPowers[0][degree-k-l]*xPowers[1][l]*xPowers[2][k];
++index;
}
}
}
return y;
}
struct Derivative
{
template<class DomainFieldType>
constexpr auto operator()(const Dune::FieldVector<DomainFieldType,3>& x) const
{
auto xPowers = std::array<std::array<RangeFieldType,maxDegree+1>,dim>{};
for(auto j : Dune::range(dim))
{
xPowers[j][0] = 1.0;
for(auto k : Dune::range(maxDegree))
xPowers[j][k+1] = xPowers[j][k]*x[j];
}
auto y = Dune::FieldMatrix<RangeFieldType,size,3>{};
std::size_t index = 0;
for(auto degree : Dune::range(maxDegree+1))
{
for(auto k : Dune::range(degree+1))
{
for (auto l : Dune::range(degree-k+1))
{
if (degree-k-l > 0)
y[index][0] = (degree-k-l)*xPowers[0][degree-k-l-1]*xPowers[1][l]*xPowers[2][k];
if (l > 0)
y[index][1] = l*xPowers[0][degree-k-l]*xPowers[1][l-1]*xPowers[2][k];
if (k > 0)
y[index][2] = k*xPowers[0][degree-k-l]*xPowers[1][l]*xPowers[2][k-1];
++index;
}
}
}
return y;
}
struct Hessian
{
template<class DomainFieldType>
constexpr auto operator()(const Dune::FieldVector<DomainFieldType,3>& x) const
{
auto xPowers = std::array<std::array<RangeFieldType,maxDegree+1>,dim>{};
for(auto j : Dune::range(dim))
Impl::computePowers<maxDegree>(x[j], xPowers[j]);
auto y = std::array<Dune::FieldMatrix<RangeFieldType,3,3>,size>{};
std::size_t index = 0;
for(auto degree : Dune::range(maxDegree+1))
{
for(auto k : Dune::range(degree+1))
{
for (auto l : Dune::range(degree-k+1))
{
// xx
if (degree-k-l-1 > 0)
y[index][0][0] = (degree-k-l)*(degree-k-l-1)*xPowers[0][degree-k-l-2]*xPowers[1][l]*xPowers[2][k];
// xy and yx
if (degree-k-l > 0 and l > 0){
y[index][0][1] = (degree-k-l)*l*xPowers[0][degree-k-l-1]*xPowers[1][l-1]*xPowers[2][k];
y[index][1][0] = y[index][0][1];
}
// yy
if (l-1 > 0)
y[index][1][1] = l*(l-1)*xPowers[0][degree-k-l]*xPowers[1][l-2]*xPowers[2][k];
// xz and zx
if (k > 0 and degree-k-l > 0)
{
y[index][0][2] = (degree-k-l)*k*xPowers[0][degree-k-l-1]*xPowers[1][l]*xPowers[2][k-1];
y[index][2][0] = y[index][0][2];
}
// yz
if (l > 0 and k > 0)
{
y[index][1][2] = l*k*xPowers[0][degree-k-l]*xPowers[1][l-1]*xPowers[2][k-1];
y[index][2][1] = y[index][1][2];
}
// zz
if (k-1 > 0)
y[index][2][2] = (k-1)*k*xPowers[0][degree-k-l]*xPowers[1][l]*xPowers[2][k-2];
++index;
}
}
}
return y;
}
};
constexpr friend auto derivative(const Derivative & d)
{
return Hessian{};
}
};
constexpr friend auto derivative(const MonomialSet& m)
{
return Derivative{};
}
};
#endif
} // namespace Dune::Functions
#endif // DUNE_FUNCTIONS_ANALYTICFUNCTIONS_MONOMIALSET_HH
Dune::Functions
Definition monomialset.hh:19
Dune::Functions::MonomialSet
Function, which evaluates all monomials up to degree maxDegree in a given coordinate.
Definition monomialset.hh:64
Dune::Functions::MonomialSet::dim
static constexpr int dim
Definition monomialset.hh:65
Dune::Functions::MonomialSet::operator()
constexpr Dune::FieldVector< RangeFieldType, size > operator()(const Dune::FieldVector< DomainFieldType, dim > &x) const
Return array of monomial evaluations.
Dune::Functions::MonomialSet::exponents
static constexpr std::array< std::array< std::size_t, dim >, size > exponents()
Return array of monomial exponents with shape size x dim.
Dune::Functions::MonomialSet::size
static constexpr int size
Definition monomialset.hh:66
Dune::Functions::MonomialSet::derivative
constexpr friend auto derivative(const MonomialSet &m)
Construct the Derivative object from a MonomialSet.
Definition monomialset.hh:145
Dune::Functions::MonomialSet::Derivative
Set of all first order derivatives of monomials up to degree maxDegree as vector of vector valued fun...
Definition monomialset.hh:99
Dune::Functions::MonomialSet::Derivative::operator()
constexpr Dune::FieldMatrix< RangeFieldType, size, dim > operator()(const Dune::FieldVector< DomainFieldType, dim > &x) const
Return array of arrays of monomial derivatives.
Dune::Functions::MonomialSet::Derivative::derivative
constexpr friend auto derivative(const Derivative &d)
Construct the Hessian object from a Derivative.
Definition monomialset.hh:135
Dune::Functions::MonomialSet::Derivative::Hessian
Set of all second order derivatives of monomials up to degree maxDegree as vector of matrix valued fu...
Definition monomialset.hh:117
Dune::Functions::MonomialSet::Derivative::Hessian::operator()
constexpr std::array< FieldMatrix< RangeFieldType, dim, dim >, size > operator()(const Dune::FieldVector< DomainFieldType, dim > &x) const
Return array of Matrices of monomial derivatives.