Select Git revision
h1functionalassembler.hh
Forked from
agnumpde / dune-fufem
1365 commits behind the upstream repository.
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
h1functionalassembler.hh 4.54 KiB
#ifndef H1_FUNCTIONAL_ASSEMBLER_HH
#define H1_FUNCTIONAL_ASSEMBLER_HH
#include <dune/common/fvector.hh>
#include <dune/common/function.hh>
#include <dune/istl/bvector.hh>
#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
#include <dune/fufem/functions/virtualgridfunction.hh>
#include <dune/fufem/assemblers/localfunctionalassembler.hh>
/** \brief Local assembler for finite element H^1-functionals
*/
template <class GridType, class T=Dune::FieldVector<double,1> >
class H1FunctionalAssembler :
public LocalFunctionalAssembler<GridType, T>
{
private:
typedef typename GridType::template Codim<0>::Geometry::GlobalCoordinate GlobalCoordinate;
typedef typename GridType::template Codim<0>::Geometry::LocalCoordinate LocalCoordinate;
static const int dim = GridType::dimension;
static const int dimworld = GridType::dimensionworld;
public:
typedef typename LocalFunctionalAssembler<GridType,T>::Element Element;
typedef typename LocalFunctionalAssembler<GridType,T>::Element::Geometry Geometry;
typedef typename LocalFunctionalAssembler<GridType,T>::LocalVector LocalVector;
typedef typename Dune::VirtualFunction<GlobalCoordinate, typename DerivativeTypefier<GlobalCoordinate,T>::DerivativeType > Function;
/** \brief constructor
*
* \param f the function representing the functional (if assembling e.g. \f$(\nabla u,\nabla v)\f$ you need to provide \f$f=\nabla u\f$)
* \param order the quadrature order (DEFAULT 2)
*/
H1FunctionalAssembler(const Function& f, int order=2) :
f_(f),
order_(order)
{}
/** \copydoc LocalFunctionalAssembler::assemble
*/
template <class TrialLocalFE>
void assemble(const Element& element, LocalVector& localVector, const TrialLocalFE& tFE)
{
typedef typename GridType::template Codim<0>::Geometry::Jacobian GeoJacobianInvTransposed;
typedef typename Function::RangeType RangeType;
typedef typename TrialLocalFE::Traits::LocalBasisType::Traits::JacobianType JacobianType;
// Make sure we got suitable shape functions
assert(tFE.type() == element.type());
// get geometry and store it
const Geometry geometry = element.geometry();
localVector = 0.0;
// get quadrature rule
// const Dune::template QuadratureRule<double, dim>& quad = Dune::template QuadratureRules<double, dim>::rule(element.type(), order_);
const Dune::template QuadratureRule<double, dim>& quad = QuadratureRuleCache<double, dim>::rule(element.type(), order_, IsRefinedLocalFiniteElement<TrialLocalFE>::value(tFE) );
// store gradients of shape functions and base functions
std::vector<JacobianType> referenceGradients(tFE.localBasis().size());
std::vector<GlobalCoordinate> gradients(tFE.localBasis().size());
// loop over quadrature points
for (size_t pt=0; pt < quad.size(); ++pt)
{
// get quadrature point
const LocalCoordinate& quadPos = quad[pt].position();
// get transposed inverse of Jacobian of transformation
const GeoJacobianInvTransposed& invJacobian = geometry.jacobianInverseTransposed(quadPos);
// get integration factor
const double integrationElement = geometry.integrationElement(quadPos);
// get gradients of shape functions
tFE.localBasis().evaluateJacobian(quadPos, referenceGradients);
// transform gradients
for (size_t i=0; i<gradients.size(); ++i)
invJacobian.mv(referenceGradients[i][0], gradients[i]);
// compute values of function
RangeType f_pos;
if (dynamic_cast<const VirtualGridFunction<GridType,RangeType>*>(&f_))
dynamic_cast<const VirtualGridFunction<GridType,RangeType>*>(&f_)->evaluateLocal(element, quadPos, f_pos);
else
f_.evaluate(geometry.global(quadPos), f_pos);
// and vector entries
for (size_t i=0; i<gradients.size(); ++i)
{
T dummy;
f_pos.mv(gradients[i],dummy);
localVector[i].axpy(quad[pt].weight()*integrationElement, dummy);
}
}
return;
}
private:
const Function& f_;
const int order_;
};
#endif