[dune-pdelab] [Patch] Assembler for l2 volume functionals
Christian Engwer
christian.engwer at uni-muenster.de
Thu Apr 10 20:55:31 CEST 2014
Hi Oli,
thanks for your patch... pushed to master
Christian
On Wed, Apr 09, 2014 at 05:55:27PM +0200, Oliver Sander wrote:
> Dear pdelab,
> please find attached a patch that implements a new local operator called L2VolumeFunctional.
> This new operator simply implements
>
> b_ = \int_\Omega f \theta_i dx
>
> for given functions f and basis vectors \theta_i. Variations of such a functional
> appear as source terms in many local operators, e.g., poisson.hh and linearelasticity.hh.
> Some code can therefore be saved by focussing on a single implementation. Also, for
> certain nonlinear energy minimization problems such l2 functionals appear without
> a corresponding matrix or differential operator.
>
> Pushing this patch to git master would be much appreciated.
>
> Thanks,
> Oliver
> From e912f77d424278706e7357d59f0844ecb8645e5b Mon Sep 17 00:00:00 2001
> From: Oliver Sander <sander at igpm.rwth-aachen.de>
> Date: Wed, 9 Apr 2014 17:45:38 +0200
> Subject: [PATCH] Add local operator that tests a given function against an FE
> space on the entire grid (aka lambda_volume)
>
> Such an operator is helpful because it appears as the source term in the
> weak form of many equations. Also, it allows to compute the integrals
> over all basis functions of an FE basis (by testing the constant 'one'
> function), which is needed for energy minimization problems involving
> superposition operators.
> ---
> dune/pdelab/localoperator/CMakeLists.txt | 3 +-
> dune/pdelab/localoperator/Makefile.am | 1 +
> dune/pdelab/localoperator/l2volumefunctional.hh | 102 ++++++++++++++++++++++++
> 3 files changed, 105 insertions(+), 1 deletion(-)
> create mode 100644 dune/pdelab/localoperator/l2volumefunctional.hh
>
> diff --git a/dune/pdelab/localoperator/CMakeLists.txt b/dune/pdelab/localoperator/CMakeLists.txt
> index 9c16088..d073a77 100644
> --- a/dune/pdelab/localoperator/CMakeLists.txt
> +++ b/dune/pdelab/localoperator/CMakeLists.txt
> @@ -18,6 +18,7 @@ set(common_HEADERS
> idefault.hh
> interface.hh
> l2.hh
> + l2volumefunctional.hh
> laplace.hh
> laplacedirichletccfv.hh
> laplacedirichletp12d.hh
> @@ -45,4 +46,4 @@ set(common_HEADERS
> # include not needed for CMake
> # include $(top_srcdir)/am/global-rules
>
> -install(FILES ${common_HEADERS} DESTINATION ${commondir})
> \ No newline at end of file
> +install(FILES ${common_HEADERS} DESTINATION ${commondir})
> diff --git a/dune/pdelab/localoperator/Makefile.am b/dune/pdelab/localoperator/Makefile.am
> index 395647e..227f850 100644
> --- a/dune/pdelab/localoperator/Makefile.am
> +++ b/dune/pdelab/localoperator/Makefile.am
> @@ -21,6 +21,7 @@ common_HEADERS = \
> idefault.hh \
> interface.hh \
> l2.hh \
> + l2volumefunctional.hh \
> laplace.hh \
> laplacedirichletccfv.hh \
> laplacedirichletp12d.hh \
> diff --git a/dune/pdelab/localoperator/l2volumefunctional.hh b/dune/pdelab/localoperator/l2volumefunctional.hh
> new file mode 100644
> index 0000000..0d2da2d
> --- /dev/null
> +++ b/dune/pdelab/localoperator/l2volumefunctional.hh
> @@ -0,0 +1,102 @@
> +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
> +// vi: set et ts=4 sw=2 sts=2:
> +#ifndef DUNE_PDELAB_LOCALOPERATOR_L2VOLUMEFUNCTIONAL_HH
> +#define DUNE_PDELAB_LOCALOPERATOR_L2VOLUMEFUNCTIONAL_HH
> +
> +#include <vector>
> +
> +#include <dune/common/exceptions.hh>
> +#include <dune/common/fvector.hh>
> +
> +#include <dune/geometry/type.hh>
> +#include <dune/geometry/quadraturerules.hh>
> +
> +#include <dune/localfunctions/common/interfaceswitch.hh>
> +
> +#include <dune/pdelab/localoperator/idefault.hh>
> +#include <dune/pdelab/localoperator/pattern.hh>
> +#include <dune/pdelab/localoperator/flags.hh>
> +
> +namespace Dune {
> + namespace PDELab {
> + //! \addtogroup LocalOperator
> + //! \ingroup PDELab
> + //! \{
> +
> + /** \brief A local operator that tests a function against a test function space defined on the entire grid
> + *
> + * It computes
> + * \f{align*}{
> + * b_i = \int_\Omega f \varphi_i\, dx
> + * \f}
> + * with conforming finite elements on all types of grids in any dimension
> + * \tparam F grid function type giving f
> + */
> + template<typename F>
> + class L2VolumeFunctional
> + : public LocalOperatorDefaultFlags
> + {
> + public:
> + // residual assembly flags
> + enum { doLambdaVolume = true };
> +
> + /** \brief Constructor
> + *
> + * \param quadOrder Order of the quadrature rule used for integrating over the element
> + */
> + L2VolumeFunctional (const F& f, unsigned int quadOrder)
> + : f_(f), quadOrder_(quadOrder)
> + {}
> +
> + // volume integral depending only on test functions
> + template<typename EG, typename LFSV, typename R>
> + void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
> + {
> + // domain and range field type
> + typedef FiniteElementInterfaceSwitch<typename LFSV::Traits::FiniteElementType> FESwitch;
> + typedef BasisInterfaceSwitch<typename FESwitch::Basis> BasisSwitch;
> +
> + typedef typename BasisSwitch::DomainField DF;
> + typedef typename BasisSwitch::RangeField RF;
> + typedef typename BasisSwitch::Range Range;
> +
> + // dimensions
> + static const int dimLocal = EG::Geometry::mydimension;
> +
> + // select quadrature rule
> + Dune::GeometryType gt = eg.geometry().type();
> + const Dune::QuadratureRule<DF,dimLocal>& rule =
> + Dune::QuadratureRules<DF,dimLocal>::rule(gt,quadOrder_);
> +
> + // loop over quadrature points
> + for (typename Dune::QuadratureRule<DF,dimLocal>::const_iterator it =
> + rule.begin(); it!=rule.end(); ++it)
> + {
> + // evaluate shape functions
> + std::vector<Range> phi(lfsv.size());
> + FESwitch::basis(lfsv.finiteElement()).
> + evaluateFunction(it->position(),phi);
> +
> + // evaluate right hand side parameter function
> + typename F::Traits::RangeType y(0.0);
> + f_.evaluate(eg.entity(),it->position(),y);
> +
> + // integrate f
> + RF factor = r.weight() * it->weight() * eg.geometry().integrationElement(it->position());
> + for (size_t i=0; i<lfsv.size(); i++)
> + r.rawAccumulate(lfsv,i,y*phi[i]*factor);
> + }
> + }
> +
> + protected:
> + const F& f_;
> +
> + // Quadrature rule order
> + unsigned int quadOrder_;
> + };
> +
> + //! \} group LocalOperator
> + } // namespace PDELab
> +} // namespace Dune
> +
> +#endif
> --
> 1.9.1
>
> _______________________________________________
> dune-pdelab mailing list
> dune-pdelab at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune-pdelab
--
Prof. Dr. Christian Engwer
Institut für Numerische und Angewandte Mathematik
Fachbereich Mathematik und Informatik der Universität Münster
Einsteinstrasse 62
48149 Münster
E-Mail christian.engwer at uni-muenster.de
Telefon +49 251 83-35067
FAX +49 251 83-32729
More information about the dune-pdelab
mailing list