Skip to content
Snippets Groups Projects
Commit 918ba8b6 authored by Carsten Gräser's avatar Carsten Gräser
Browse files

Reverse dependency of basisgridfunction.hh and vintagebasisgridfunction.hh 

This moves the `Local(Derivative)Evaluator` implementations from
`basisgridfunction.hh` to `vintagebasisgridfunction.hh`. Hence
the former now includes the latter. This was the other way around before.

By this change one can include `vintagebasisgridfunction.hh` without
seeing deprecation warnings because `dune/common/function.hh` is implicitly
included. Furthermore it allows to remove `basisgridfunction.hh` before
`vintagebasisgridfunction.hh`.
parent 9514ec94
No related branches found
No related tags found
1 merge request!146Stop using Dune::VirtualFunction interface
Hide whitespace changes
Inline Side-by-side
dune/fufem/functions/basisgridfunction.hh
+ 3
− 120
View file @ 918ba8b6
......@@ -11,11 +11,9 @@
#include <dune/istl/bvector.hh>
#include <dune/matrix-vector/axpy.hh>
#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
#include <dune/fufem/functions/virtualgridfunction.hh>
#include <dune/fufem/functions/vintagebasisgridfunction.hh>
template<class B, class C> class BasisGridFunction;
......@@ -25,121 +23,6 @@ template<class B, class C> class BasisGridFunction;
template<class Basis, class Coefficients, class RangeType>
struct LocalEvaluator
{
template<class Element, class LocalDomainType>
static void evaluateLocal(const Basis& basis,
const Coefficients& coefficients,
const Element& e,
const LocalDomainType& x,
RangeType& y)
{
typedef typename Basis::LinearCombination LinearCombination;
typedef typename Basis::LocalFiniteElement::Traits::LocalBasisType LB;
typedef typename LB::Traits::RangeType LBRangeType;
const LB& localBasis = basis.getLocalFiniteElement(e).localBasis();
std::vector<LBRangeType> values;
localBasis.evaluateFunction(x, values);
y = 0;
for (size_t i=0; i<localBasis.size(); ++i)
{
int index = basis.index(e, i);
if (basis.isConstrained(index))
{
const LinearCombination& constraints = basis.constraints(index);
for (size_t w=0; w<constraints.size(); ++w)
Dune::MatrixVector::addProduct(y, values[i] * constraints[w].factor, coefficients[constraints[w].index]);
}
else
Dune::MatrixVector::addProduct(y, values[i], coefficients[index]);
}
}
};
template<class Basis, class Coefficients, class RangeType>
struct LocalDerivativeEvaluator
{
template<class Element, class LocalDomainType, class DerivativeType>
static void evaluateDerivativeLocal(const Basis& basis,
const Coefficients& c,
const Element& e,
const LocalDomainType& x,
DerivativeType& d)
{
DUNE_THROW(Dune::NotImplemented, "evaluateDerivativeLocal not implemented for this RangeType");
}
};
template<class Basis, class Coefficients, class RangeFieldType, int m>
struct LocalDerivativeEvaluator<Basis, Coefficients, typename Dune::FieldVector<RangeFieldType,m> >
{
//! Compute y+= a*A*B
template<class Weight, class Jacobian, class DerivativeType>
static void gaxpy(const Weight& a,
const typename Coefficients::value_type& A,
const Jacobian& B,
DerivativeType& y)
{
for (int i=0; i<DerivativeType::rows; ++i)
for (int j=0; j<DerivativeType::cols; ++j)
y[i][j] += a * A[i] * B[0][j];
}
//! Compute y+= A*B
template<class Jacobian, class DerivativeType>
static void gaxpy(const typename Coefficients::value_type& A,
const Jacobian& B,
DerivativeType& y)
{
for (int i=0; i<DerivativeType::rows; ++i)
for (int j=0; j<DerivativeType::cols; ++j)
y[i][j] += A[i] * B[0][j];
}
template<class Element, class LocalDomainType, class DerivativeType>
static void evaluateDerivativeLocal(const Basis& basis,
const Coefficients& coefficients,
const Element& e,
const LocalDomainType& x,
DerivativeType& d)
{
typedef typename Basis::LinearCombination LinearCombination;
typedef typename Basis::LocalFiniteElement::Traits::LocalBasisType LB;
typedef typename LB::Traits::JacobianType LBJacobianType;
const LB& localBasis = basis.getLocalFiniteElement(e).localBasis();
std::vector<LBJacobianType> localJacobians;
localBasis.evaluateJacobian(x, localJacobians);
DerivativeType localD(0.0);
for (size_t i=0; i<localBasis.size(); ++i)
{
int index = basis.index(e, i);
if (basis.isConstrained(index))
{
const LinearCombination& constraints = basis.constraints(index);
for (size_t w=0; w<constraints.size(); ++w)
gaxpy(constraints[w].factor, coefficients[constraints[w].index], localJacobians[i], localD);
}
else
gaxpy(coefficients[index], localJacobians[i], localD);
}
typename Element::Geometry::JacobianInverseTransposed J = e.geometry().jacobianInverseTransposed(x);
for (int j=0; j<DerivativeType::rows; ++j)
J.mv(localD[j], d[j]);
}
};
/**
......@@ -229,7 +112,7 @@ class BasisGridFunction :
*/
virtual void evaluateLocal(const Element& e, const LocalDomainType& x, RangeType& y) const
{
LocalEvaluator<B, CoefficientVector, RangeType>::evaluateLocal(basis_, coefficients_, e, x, y);
Dune::Fufem::Impl::LocalEvaluator<B, CoefficientVector, RangeType>::evaluateLocal(basis_, coefficients_, e, x, y);
}
......@@ -239,7 +122,7 @@ class BasisGridFunction :
*/
void evaluateDerivativeLocal(const Element& e, const LocalDomainType& x, DerivativeType& d) const
{
LocalDerivativeEvaluator<B, CoefficientVector, RangeType>::evaluateDerivativeLocal(basis_, coefficients_, e, x, d);
Dune::Fufem::Impl::LocalDerivativeEvaluator<B, CoefficientVector, DerivativeType>::evaluateLocal(basis_, coefficients_, e, x, d);
}
......
dune/fufem/functions/vintagebasisgridfunction.hh
+ 127
− 2
View file @ 918ba8b6
......@@ -10,15 +10,140 @@
#include <dune/common/shared_ptr.hh>
#include <dune/common/exceptions.hh>
#include <dune/matrix-vector/axpy.hh>
#include <dune/functions/gridfunctions/gridviewentityset.hh>
#include <dune/functions/gridfunctions/gridfunction.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
namespace Dune::Fufem {
namespace Impl {
template<class Basis, class Coefficients, class DerivativeRangeType>
struct LocalDerivativeEvaluator;
template<class Basis, class Coefficients, class RangeType>
struct LocalEvaluator
{
template<class Element, class LocalDomainType>
static void evaluateLocal(const Basis& basis,
const Coefficients& coefficients,
const Element& e,
const LocalDomainType& x,
RangeType& y)
{
typedef typename Basis::LinearCombination LinearCombination;
typedef typename Basis::LocalFiniteElement::Traits::LocalBasisType LB;
typedef typename LB::Traits::RangeType LBRangeType;
const LB& localBasis = basis.getLocalFiniteElement(e).localBasis();
std::vector<LBRangeType> values;
localBasis.evaluateFunction(x, values);
y = 0;
for (size_t i=0; i<localBasis.size(); ++i)
{
int index = basis.index(e, i);
if (basis.isConstrained(index))
{
const LinearCombination& constraints = basis.constraints(index);
for (size_t w=0; w<constraints.size(); ++w)
Dune::MatrixVector::addProduct(y, values[i] * constraints[w].factor, coefficients[constraints[w].index]);
}
else
Dune::MatrixVector::addProduct(y, values[i], coefficients[index]);
}
}
};
template<class Basis, class Coefficients, class DerivativeRangeType>
struct LocalDerivativeEvaluator
{
template<class Element, class LocalDomainType>
static void evaluateLocal(const Basis& basis,
const Coefficients& c,
const Element& e,
const LocalDomainType& x,
DerivativeRangeType& d)
{
DUNE_THROW(Dune::NotImplemented, "LocalDerivativeEvaluator not implemented for this DerivativeRangeType");
}
};
template<class Basis, class Coefficients, class RangeFieldType, int n, int m>
struct LocalDerivativeEvaluator<Basis, Coefficients, typename Dune::FieldMatrix<RangeFieldType,n,m> >
{
using DerivativeRangeType = typename Dune::FieldMatrix<RangeFieldType,n,m>;
//! Compute y+= a*A*B
template<class Weight, class Jacobian>
static void gaxpy(const Weight& a,
const typename Coefficients::value_type& A,
const Jacobian& B,
DerivativeRangeType& y)
{
for (int i=0; i<DerivativeRangeType::rows; ++i)
for (int j=0; j<DerivativeRangeType::cols; ++j)
y[i][j] += a * A[i] * B[0][j];
}
//! Compute y+= A*B
template<class Jacobian>
static void gaxpy(const typename Coefficients::value_type& A,
const Jacobian& B,
DerivativeRangeType& y)
{
for (int i=0; i<DerivativeRangeType::rows; ++i)
for (int j=0; j<DerivativeRangeType::cols; ++j)
y[i][j] += A[i] * B[0][j];
}
template<class Element, class LocalDomainType>
static void evaluateLocal(const Basis& basis,
const Coefficients& coefficients,
const Element& e,
const LocalDomainType& x,
DerivativeRangeType& d)
{
typedef typename Basis::LinearCombination LinearCombination;
typedef typename Basis::LocalFiniteElement::Traits::LocalBasisType LB;
typedef typename LB::Traits::JacobianType LBJacobianType;
const LB& localBasis = basis.getLocalFiniteElement(e).localBasis();
std::vector<LBJacobianType> localJacobians;
localBasis.evaluateJacobian(x, localJacobians);
DerivativeRangeType localD(0.0);
for (size_t i=0; i<localBasis.size(); ++i)
{
int index = basis.index(e, i);
if (basis.isConstrained(index))
{
const LinearCombination& constraints = basis.constraints(index);
for (size_t w=0; w<constraints.size(); ++w)
gaxpy(constraints[w].factor, coefficients[constraints[w].index], localJacobians[i], localD);
}
else
gaxpy(coefficients[index], localJacobians[i], localD);
}
typename Element::Geometry::JacobianInverseTransposed J = e.geometry().jacobianInverseTransposed(x);
for (int j=0; j<DerivativeRangeType::rows; ++j)
J.mv(localD[j], d[j]);
}
};
/**
* \brief A dune-functions grid function implemented using a dune-fufem basis
*
......@@ -225,7 +350,7 @@ auto makeVintageBasisGridFunction(B&& basis, V&& v)
using Basis = std::decay_t<B>;
using Vector = std::decay_t<V>;
auto localEvaluator = [](const auto& b, const auto& v, const auto&e, const auto& x, auto& y) {
LocalEvaluator<Basis, Vector, R>::evaluateLocal(b, v, e, x, y);
Dune::Fufem::Impl::LocalEvaluator<Basis, Vector, R>::evaluateLocal(b, v, e, x, y);
};
using LocalEvaluator = decltype(localEvaluator);
return Impl::VintageBasisGridFunction<Basis, Vector, LocalEvaluator, R>(std::forward<B>(basis), std::forward<V>(v), std::move(localEvaluator));
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment