Loading...
Searching...
No Matches
globalvaluedlocalfiniteelement.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_GLOBALVALUEDLOCALFINITEELEMENT_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_GLOBALVALUEDLOCALFINITEELEMENT_HH
#include <array>
#include <numeric>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/math.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/localfunctions/common/localbasis.hh>
#include <dune/localfunctions/common/localfiniteelementtraits.hh>
#include <dune/functions/common/densevectorview.hh>
namespace Dune::Functions::Impl
{
struct ContravariantPiolaTransformator
{
template<typename Values, typename LocalCoordinate, typename Geometry>
static auto apply(Values& values,
const LocalCoordinate& xi,
const Geometry& geometry)
{
auto jacobianTransposed = geometry.jacobianTransposed(xi);
auto integrationElement = geometry.integrationElement(xi);
for (auto& value : values)
{
auto tmp = value;
// Here we cannot directly use
// jacobianTransposed.mtv(tmp, value);
// because mtv expects the DenseVector interface.
auto tmpDenseVector = Impl::DenseVectorView(tmp);
auto valueDenseVector = Impl::DenseVectorView(value);
jacobianTransposed.mtv(tmpDenseVector, valueDenseVector);
value /= integrationElement;
}
}
template<typename Gradients, typename LocalCoordinate, typename Geometry>
static auto applyJacobian(Gradients& gradients,
const LocalCoordinate& xi,
const Geometry& geometry)
{
auto jacobianTransposed = geometry.jacobianTransposed(xi);
auto integrationElement = geometry.integrationElement(xi);
for (auto& gradient : gradients)
{
auto tmp = gradient;
gradient = 0;
for (size_t k=0; k<gradient.M(); k++)
for (size_t l=0; l<tmp.N(); l++)
// Use sparseRange because jacobianTransposed may be a sparse DiagonalMatrix
for(auto&& [jacobianTransposed_l_j, j] : sparseRange(jacobianTransposed[l]))
gradient[j][k] += jacobianTransposed_l_j * tmp[l][k];
gradient /= integrationElement;
}
}
template<class Function, class LocalCoordinate, class Element>
class LocalValuedFunction
{
const Function& f_;
const Element& element_;
using LocalValue = LocalCoordinate;
public:
LocalValuedFunction(const Function& f, const Element& element)
: f_(f), element_(element)
{}
auto operator()(const LocalCoordinate& xi) const
{
auto globalValue = f_(xi);
// Apply the inverse Piola transform
auto jacobianInverseTransposed = element_.geometry().jacobianInverseTransposed(xi);
auto integrationElement = element_.geometry().integrationElement(xi);
auto localValue = LocalValue{};
// Here we cannot directly use
// jacobianInverseTransposed.mtv(globalValue, localValue);
// because mtv expects the DenseVector interface.
auto globalValueDenseVector = Impl::DenseVectorView(globalValue);
auto localValueDenseVector = Impl::DenseVectorView(localValue);
jacobianInverseTransposed.mtv(globalValueDenseVector, localValueDenseVector);
localValue *= integrationElement;
return localValue;
}
};
};
struct CovariantPiolaTransformator
{
template<typename Values, typename LocalCoordinate, typename Geometry>
static auto apply(Values& values,
const LocalCoordinate& xi,
const Geometry& geometry)
{
auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(xi);
for (auto& value : values)
{
auto tmp = value;
// Here we cannot directly use
// jacobianInverseTransposed.mv(tmp, value);
// because mv expects the DenseVector interface.
auto tmpDenseVector = Impl::DenseVectorView(tmp);
auto valueDenseVector = Impl::DenseVectorView(value);
jacobianInverseTransposed.mv(tmpDenseVector, valueDenseVector);
}
}
template<typename Gradients, typename LocalCoordinate, typename Geometry>
static auto applyJacobian(Gradients& gradients,
const LocalCoordinate& xi,
const Geometry& geometry)
{
auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(xi);
for (auto& gradient : gradients)
{
auto tmp = gradient;
gradient = 0;
for (size_t j=0; j<gradient.N(); j++)
for (size_t k=0; k<gradient.M(); k++)
// Use sparseRange because jacobianTransposed may be a sparse DiagonalMatrix
for(auto&& [jacobianInverseTransposed_j_l, l] : sparseRange(jacobianInverseTransposed[j]))
gradient[j][k] += jacobianInverseTransposed_j_l * tmp[l][k];
}
}
template<class Function, class LocalCoordinate, class Element>
class LocalValuedFunction
{
const Function& f_;
const Element& element_;
public:
LocalValuedFunction(const Function& f, const Element& element)
: f_(f), element_(element)
{}
auto operator()(const LocalCoordinate& xi) const
{
auto globalValue = f_(xi);
// Apply the inverse Piola transform
auto jacobianTransposed = element_.geometry().jacobianTransposed(xi);
auto localValue = globalValue;
// Here we cannot directly use
// jacobianTransposed.mv(globalValue, localValue);
// because mv expects the DenseVector interface.
auto localValueDenseVector = Impl::DenseVectorView(localValue);
auto globalValueDenseVector = Impl::DenseVectorView(globalValue);
jacobianTransposed.mv(globalValueDenseVector, localValueDenseVector);
return localValue;
}
};
};
template<class Transformator, class LocalValuedLocalBasis, class Element>
class GlobalValuedLocalBasis
{
public:
using Traits = typename LocalValuedLocalBasis::Traits;
void bind(const LocalValuedLocalBasis& localValuedLocalBasis, const Element& element)
{
localValuedLocalBasis_ = &localValuedLocalBasis;
element_ = &element;
}
auto size() const
{
return localValuedLocalBasis_->size();
}
void evaluateFunction(const typename Traits::DomainType& x,
std::vector<typename Traits::RangeType>& out) const
{
localValuedLocalBasis_->evaluateFunction(x,out);
Transformator::apply(out, x, element_->geometry());
}
void evaluateJacobian(const typename Traits::DomainType& x,
std::vector<typename Traits::JacobianType>& out) const
{
localValuedLocalBasis_->evaluateJacobian(x,out);
Transformator::applyJacobian(out, x, element_->geometry());
}
void partial(const std::array<unsigned int,2>& order,
const typename Traits::DomainType& x,
std::vector<typename Traits::RangeType>& out) const
{
auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
if (totalOrder == 0) {
evaluateFunction(x, out);
} else if (totalOrder == 1) {
auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
out.resize(size());
// TODO: The following is wasteful: We compute the full Jacobian and then return
// only a part of it. While we need the full Jacobian of the underlying local-valued LFE,
// it should be possible to compute only a partial Piola transform for the requested
// partial derivatives.
std::vector<typename Traits::JacobianType> fullJacobian;
localValuedLocalBasis_->evaluateJacobian(x,fullJacobian);
Transformator::applyJacobian(fullJacobian, x, element_->geometry());
for (std::size_t i=0; i<out.size(); i++)
for (std::size_t j=0; j<out[i].size(); j++)
out[i][j] = fullJacobian[i][j][direction];
} else
DUNE_THROW(NotImplemented, "Partial derivatives of order 2 or higher");
}
auto order() const
{
return localValuedLocalBasis_->order();
}
const LocalValuedLocalBasis* localValuedLocalBasis_;
const Element* element_;
};
template<class Transformator, class LocalValuedLocalInterpolation, class Element>
class GlobalValuedLocalInterpolation
{
public:
void bind(const LocalValuedLocalInterpolation& localValuedLocalInterpolation, const Element& element)
{
localValuedLocalInterpolation_ = &localValuedLocalInterpolation;
element_ = &element;
}
template<typename F, typename C>
void interpolate (const F& f, std::vector<C>& out) const
{
using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
typename Transformator::template LocalValuedFunction<F,LocalCoordinate,Element> localValuedFunction(f, *element_);
localValuedLocalInterpolation_->interpolate(localValuedFunction, out);
}
private:
const LocalValuedLocalInterpolation* localValuedLocalInterpolation_;
const Element* element_;
};
template<class Transformator, class LocalValuedLFE, class Element>
class GlobalValuedLocalFiniteElement
{
using LocalBasis = GlobalValuedLocalBasis<Transformator,
typename LocalValuedLFE::Traits::LocalBasisType,
Element>;
using LocalInterpolation = GlobalValuedLocalInterpolation<Transformator,
typename LocalValuedLFE::Traits::LocalInterpolationType,
Element>;
public:
using Traits = LocalFiniteElementTraits<LocalBasis,
typename LocalValuedLFE::Traits::LocalCoefficientsType,
LocalInterpolation>;
GlobalValuedLocalFiniteElement() {}
void bind(const LocalValuedLFE& localValuedLFE, const Element& element)
{
globalValuedLocalBasis_.bind(localValuedLFE.localBasis(), element);
globalValuedLocalInterpolation_.bind(localValuedLFE.localInterpolation(), element);
localValuedLFE_ = &localValuedLFE;
}
const typename Traits::LocalBasisType& localBasis() const
{
return globalValuedLocalBasis_;
}
const typename Traits::LocalCoefficientsType& localCoefficients() const
{
return localValuedLFE_->localCoefficients();
}
const typename Traits::LocalInterpolationType& localInterpolation() const
{
return globalValuedLocalInterpolation_;
}
std::size_t size() const
{
return localValuedLFE_->size();
}
GeometryType type() const
{
return localValuedLFE_->type();
}
private:
typename Traits::LocalBasisType globalValuedLocalBasis_;
typename Traits::LocalInterpolationType globalValuedLocalInterpolation_;
const LocalValuedLFE* localValuedLFE_;
};
} // namespace Dune::Functions::Impl
#endif // DUNE_FUNCTIONS_FUNCTIONSPACEBASES_GLOBALVALUEDLOCALFINITEELEMENT_HH
densevectorview.hh