[dune-pdelab] Patches for local operators
Christian Engwer
christian.engwer at uni-muenster.de
Sat Apr 5 00:05:05 CEST 2014
Hi Oli,
thanks for the patches.
In your updats to the existing operators you instantiate new instances
of Laplace in every call to alpha_volume / jacobian_volume.
I'd actually prefer to have the original Laplace, either as a base
class, or as a private member.
I added the laplace operator and hope you provide updates for the
rewrite of existing operators.
Ciao
Christian
On Sat, Mar 29, 2014 at 12:11:02PM +0100, Oliver Sander wrote:
> Dear pdelab,
> please find attached five patches for dune-pdelab
>
> The first adds a new LocalOperator "Laplace" which does nothing except assemble
> the Laplace matrix. It does not know about boundary values, volume terms or
> constraints, and caters to people who value simplicity.
>
> The next two patches use the new operator to simplify the implementation of the
> Poisson operator
>
> The last two patches are unrelated; they do some minor cleanup of the LinearElasticity
> operator.
>
> I'd be happy if these patches could be applied to the dune-pdelab master branch.
>
> Thanks,
> Oliver
> From 6f83d835cf60c4e8286067d8f4fdb36e8b61c537 Mon Sep 17 00:00:00 2001
> From: Oliver Sander <sander at igpm.rwth-aachen.de>
> Date: Sat, 29 Mar 2014 10:38:26 +0100
> Subject: [PATCH 1/5] Add a local operator for the Laplace problem
>
> This operator only assembles the Laplace stiffness matrix, and nothing
> else. No volume terms, no Neumann terms, no constraints.
>
> Supported is:
> - explicitly assembling the matrix (aka "jacobian_volume")
> - matrix-vectors products without explicit matrix assembly
> (aka "alpha_volume")
> ---
> dune/pdelab/localoperator/CMakeLists.txt | 1 +
> dune/pdelab/localoperator/Makefile.am | 1 +
> dune/pdelab/localoperator/laplace.hh | 178 +++++++++++++++++++++++++++++++
> 3 files changed, 180 insertions(+)
> create mode 100644 dune/pdelab/localoperator/laplace.hh
>
> diff --git a/dune/pdelab/localoperator/CMakeLists.txt b/dune/pdelab/localoperator/CMakeLists.txt
> index 8c9a6d9..9c16088 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
> + laplace.hh
> laplacedirichletccfv.hh
> laplacedirichletp12d.hh
> linearelasticity.hh
> diff --git a/dune/pdelab/localoperator/Makefile.am b/dune/pdelab/localoperator/Makefile.am
> index c4aac27..395647e 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 \
> + laplace.hh \
> laplacedirichletccfv.hh \
> laplacedirichletp12d.hh \
> linearelasticity.hh \
> diff --git a/dune/pdelab/localoperator/laplace.hh b/dune/pdelab/localoperator/laplace.hh
> new file mode 100644
> index 0000000..68ee8c9
> --- /dev/null
> +++ b/dune/pdelab/localoperator/laplace.hh
> @@ -0,0 +1,178 @@
> +// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
> +// vi: set et ts=4 sw=2 sts=2:
> +
> +#ifndef DUNE_PDELAB_LOCALOPERATOR_LAPLACE_HH
> +#define DUNE_PDELAB_LOCALOPERATOR_LAPLACE_HH
> +
> +#include <vector>
> +
> +#include <dune/common/fvector.hh>
> +#include <dune/common/static_assert.hh>
> +
> +#include <dune/geometry/type.hh>
> +#include <dune/geometry/quadraturerules.hh>
> +
> +#include <dune/localfunctions/common/interfaceswitch.hh>
> +
> +#include <dune/pdelab/localoperator/pattern.hh>
> +#include <dune/pdelab/localoperator/flags.hh>
> +
> +namespace Dune {
> + namespace PDELab {
> + //! \addtogroup LocalOperator
> + //! \ingroup PDELab
> + //! \{
> +
> + /** a local operator for solving the Laplace equation
> + *
> + * \f{align*}{
> + * - \Delta u &=& 0 \mbox{ in } \Omega, \\
> + * -\nabla u \cdot \nu &=& 0 \mbox{ on } \partial\Omega_N \\
> + * \f}
> + * with conforming finite elements on all types of grids in any dimension.
> + *
> + * In other words, it only assembles the Laplace matrix.
> + *
> + */
> + class Laplace
> + : public FullVolumePattern,
> + public LocalOperatorDefaultFlags
> + {
> + public:
> + // pattern assembly flags
> + enum { doPatternVolume = true };
> +
> + // residual assembly flags
> + enum { doAlphaVolume = true };
> +
> + /** \brief Constructor
> + *
> + * \param quadOrder Order of the quadrature rule used for integrating over the element
> + */
> + Laplace (unsigned int quadOrder)
> + : quadOrder_(quadOrder)
> + {}
> +
> + /** \brief Compute Laplace matrix times a given vector for one element
> + *
> + * This is used for matrix-free algorithms for the Laplace equation
> + *
> + * \param [in] eg The grid element we are assembling on
> + * \param [in] lfsu Local ansatz function space basis
> + * \param [in] lfsv Local test function space basis
> + * \param [in] x Input vector
> + * \param [out] r The product of the Laplace matrix times x
> + */
> + template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
> + void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
> + {
> + // domain and range field type
> + typedef FiniteElementInterfaceSwitch<
> + typename LFSU::Traits::FiniteElementType
> + > FESwitch;
> + typedef BasisInterfaceSwitch<
> + typename FESwitch::Basis
> + > BasisSwitch;
> + typedef typename BasisSwitch::DomainField DF;
> + typedef typename BasisSwitch::RangeField RF;
> +
> + // dimensions
> + static const int dimLocal = EG::Geometry::mydimension;
> + static const int dimGlobal = EG::Geometry::coorddimension;
> +
> + // 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 gradient of shape functions
> + // (we assume Galerkin method lfsu=lfsv)
> + std::vector<Dune::FieldMatrix<RF,1,dimGlobal> >
> + gradphiu(lfsu.size());
> + BasisSwitch::gradient(FESwitch::basis(lfsu.finiteElement()),
> + eg.geometry(), it->position(), gradphiu);
> + std::vector<Dune::FieldMatrix<RF,1,dimGlobal> >
> + gradphiv(lfsv.size());
> + BasisSwitch::gradient(FESwitch::basis(lfsv.finiteElement()),
> + eg.geometry(), it->position(), gradphiv);
> +
> + // compute gradient of u
> + Dune::FieldVector<RF,dimGlobal> gradu(0.0);
> + for (size_t i=0; i<lfsu.size(); i++)
> + gradu.axpy(x(lfsu,i),gradphiu[i][0]);
> +
> + // integrate grad u * grad phi_i
> + RF factor = r.weight() * it->weight() * eg.geometry().integrationElement(it->position());
> + for (size_t i=0; i<lfsv.size(); i++)
> + r.rawAccumulate(lfsv,i,(gradu*gradphiv[i][0])*factor);
> + }
> + }
> +
> + /** \brief Compute the Laplace stiffness matrix for the element given in 'eg'
> + *
> + * \tparam M Type of the element stiffness matrix
> + *
> + * \param [in] eg The grid element we are assembling on
> + * \param [in] lfsu Local ansatz function space basis
> + * \param [in] lfsv Local test function space basis
> + * \param [in] x Current configuration; gets ignored for linear problems like this one
> + * \param [out] matrix Element stiffness matrix
> + */
> + template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
> + void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, M & matrix) const
> + {
> + // Switches between local and global interface
> + typedef FiniteElementInterfaceSwitch<
> + typename LFSU::Traits::FiniteElementType
> + > FESwitch;
> + typedef BasisInterfaceSwitch<
> + typename FESwitch::Basis
> + > BasisSwitch;
> +
> + // domain and range field type
> + typedef typename BasisSwitch::DomainField DF;
> + typedef typename BasisSwitch::RangeField RF;
> + typedef typename LFSU::Traits::SizeType size_type;
> +
> + // dimensions
> + const int dim = EG::Geometry::dimension;
> +
> + // select quadrature rule
> + Dune::GeometryType gt = eg.geometry().type();
> + const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,quadOrder_);
> +
> + // loop over quadrature points
> + for (typename QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
> + {
> + std::vector<Dune::FieldMatrix<RF,1,dim> > gradphi(lfsu.size());
> + BasisSwitch::gradient(FESwitch::basis(lfsu.finiteElement()),
> + eg.geometry(), it->position(), gradphi);
> +
> + // geometric weight
> + RF factor = it->weight() * eg.geometry().integrationElement(it->position());
> +
> + for (size_type i=0; i<lfsu.size(); i++)
> + {
> + for (size_type j=0; j<lfsv.size(); j++)
> + {
> + // integrate grad u * grad phi
> + matrix.accumulate(lfsv,j,lfsu,i, gradphi[i][0] * gradphi[j][0] * factor);
> + }
> + }
> + }
> + }
> +
> + protected:
> + // Quadrature rule order
> + unsigned int quadOrder_;
> + };
> +
> + //! \} group LocalOperator
> + } // namespace PDELab
> +} // namespace Dune
> +
> +#endif
> --
> 1.9.0
>
> From d1dba9cf7dae17c7e09c7721ad1b1e0e2b96ad2a Mon Sep 17 00:00:00 2001
> From: Oliver Sander <sander at igpm.rwth-aachen.de>
> Date: Sat, 29 Mar 2014 10:50:31 +0100
> Subject: [PATCH 2/5] Use the new Laplace assembler (in laplace.hh) to
> implement alpha_volume method
>
> ---
> dune/pdelab/localoperator/poisson.hh | 48 +++---------------------------------
> 1 file changed, 4 insertions(+), 44 deletions(-)
>
> diff --git a/dune/pdelab/localoperator/poisson.hh b/dune/pdelab/localoperator/poisson.hh
> index a3c2122..4eef0af 100644
> --- a/dune/pdelab/localoperator/poisson.hh
> +++ b/dune/pdelab/localoperator/poisson.hh
> @@ -15,6 +15,8 @@
>
> #include <dune/localfunctions/common/interfaceswitch.hh>
>
> +#include <dune/pdelab/localoperator/laplace.hh>
> +
> #include"defaultimp.hh"
> #include"idefault.hh"
> #include"pattern.hh"
> @@ -66,50 +68,8 @@ namespace Dune {
> template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
> void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
> {
> - // domain and range field type
> - typedef FiniteElementInterfaceSwitch<
> - typename LFSU::Traits::FiniteElementType
> - > FESwitch;
> - typedef BasisInterfaceSwitch<
> - typename FESwitch::Basis
> - > BasisSwitch;
> - typedef typename BasisSwitch::DomainField DF;
> - typedef typename BasisSwitch::RangeField RF;
> -
> - // dimensions
> - static const int dimLocal = EG::Geometry::mydimension;
> - static const int dimGlobal = EG::Geometry::coorddimension;
> -
> - // 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 gradient of shape functions
> - // (we assume Galerkin method lfsu=lfsv)
> - std::vector<Dune::FieldMatrix<RF,1,dimGlobal> >
> - gradphiu(lfsu.size());
> - BasisSwitch::gradient(FESwitch::basis(lfsu.finiteElement()),
> - eg.geometry(), it->position(), gradphiu);
> - std::vector<Dune::FieldMatrix<RF,1,dimGlobal> >
> - gradphiv(lfsv.size());
> - BasisSwitch::gradient(FESwitch::basis(lfsv.finiteElement()),
> - eg.geometry(), it->position(), gradphiv);
> -
> - // compute gradient of u
> - Dune::FieldVector<RF,dimGlobal> gradu(0.0);
> - for (size_t i=0; i<lfsu.size(); i++)
> - gradu.axpy(x(lfsu,i),gradphiu[i][0]);
> -
> - // integrate grad u * grad phi_i
> - RF factor = r.weight() * it->weight() * eg.geometry().integrationElement(it->position());
> - for (size_t i=0; i<lfsv.size(); i++)
> - r.rawAccumulate(lfsv,i,(gradu*gradphiv[i][0])*factor);
> - }
> + Laplace laplace(quadOrder_);
> + laplace.alpha_volume(eg, lfsu, x, lfsv, r);
> }
>
> // volume integral depending only on test functions
> --
> 1.9.0
>
> From 6f93b136f04d66f75b63e31750ab89941001fd67 Mon Sep 17 00:00:00 2001
> From: Oliver Sander <sander at igpm.rwth-aachen.de>
> Date: Sat, 29 Mar 2014 11:24:12 +0100
> Subject: [PATCH 3/5] Implement the matrix assemble (i.e., jacobian_volume),
> instead of relying on an FD approximation
>
> This is very short: we simply use the implementation in laplace.hh
> ---
> dune/pdelab/localoperator/poisson.hh | 17 ++++++++++++++++-
> 1 file changed, 16 insertions(+), 1 deletion(-)
>
> diff --git a/dune/pdelab/localoperator/poisson.hh b/dune/pdelab/localoperator/poisson.hh
> index 4eef0af..a6bd974 100644
> --- a/dune/pdelab/localoperator/poisson.hh
> +++ b/dune/pdelab/localoperator/poisson.hh
> @@ -42,7 +42,6 @@ namespace Dune {
> */
> template<typename F, typename B, typename J>
> class Poisson : public NumericalJacobianApplyVolume<Poisson<F,B,J> >,
> - public NumericalJacobianVolume<Poisson<F,B,J> >,
> public FullVolumePattern,
> public LocalOperatorDefaultFlags
> {
> @@ -72,6 +71,22 @@ namespace Dune {
> laplace.alpha_volume(eg, lfsu, x, lfsv, r);
> }
>
> + /** \brief Compute the Laplace stiffness matrix for the element given in 'eg'
> + *
> + * \tparam M Type of the element stiffness matrix
> + *
> + * \param [in] eg The grid element we are assembling on
> + * \param [in] lfsu Local ansatz function space basis
> + * \param [in] lfsv Local test function space basis
> + * \param [in] x Current configuration; gets ignored for linear problems like this one
> + * \param [out] matrix Element stiffness matrix
> + */
> + template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
> + void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, M & matrix) const
> + {
> + Laplace laplace(quadOrder_);
> + laplace.jacobian_volume(eg, lfsu, x, lfsv, matrix);
> + }
> // 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
> --
> 1.9.0
>
> From c9d318061c27725294bdf8b4bb53a8517ec60a22 Mon Sep 17 00:00:00 2001
> From: Oliver Sander <sander at igpm.rwth-aachen.de>
> Date: Sat, 29 Mar 2014 11:40:20 +0100
> Subject: [PATCH 4/5] [cleanup] Do not evaluate shape function values. They
> are never used.
>
> ---
> dune/pdelab/localoperator/linearelasticity.hh | 12 ------------
> 1 file changed, 12 deletions(-)
>
> diff --git a/dune/pdelab/localoperator/linearelasticity.hh b/dune/pdelab/localoperator/linearelasticity.hh
> index b0e5ee8..8e001f1 100644
> --- a/dune/pdelab/localoperator/linearelasticity.hh
> +++ b/dune/pdelab/localoperator/linearelasticity.hh
> @@ -70,8 +70,6 @@ namespace Dune {
> Traits::LocalBasisType::Traits::RangeFieldType RF;
> typedef typename LFSU_SUB::Traits::FiniteElementType::
> Traits::LocalBasisType::Traits::JacobianType JacobianType;
> - typedef typename LFSU_SUB::Traits::FiniteElementType::
> - Traits::LocalBasisType::Traits::RangeType RangeType;
>
> typedef typename LFSU_SUB::Traits::SizeType size_type;
>
> @@ -87,10 +85,6 @@ namespace Dune {
> // loop over quadrature points
> for (typename QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
> {
> - // evaluate basis functions
> - std::vector<RangeType> phi(lfsu.child(0).size());
> - lfsu.child(0).finiteElement().localBasis().evaluateFunction(it->position(),phi);
> -
> // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
> std::vector<JacobianType> js(lfsu.child(0).size());
> lfsu.child(0).finiteElement().localBasis().evaluateJacobian(it->position(),js);
> @@ -154,8 +148,6 @@ namespace Dune {
> Traits::LocalBasisType::Traits::RangeFieldType RF;
> typedef typename LFSU::Traits::FiniteElementType::
> Traits::LocalBasisType::Traits::JacobianType JacobianType;
> - typedef typename LFSU::Traits::FiniteElementType::
> - Traits::LocalBasisType::Traits::RangeType RangeType;
>
> typedef typename LFSU::Traits::SizeType size_type;
>
> @@ -171,10 +163,6 @@ namespace Dune {
> // loop over quadrature points
> for (typename QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
> {
> - // evaluate basis functions
> - std::vector<RangeType> phi(lfsu_hat.child(0).size());
> - lfsu_hat.child(0).finiteElement().localBasis().evaluateFunction(it->position(),phi);
> -
> // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
> std::vector<JacobianType> js(lfsu_hat.child(0).size());
> lfsu_hat.child(0).finiteElement().localBasis().evaluateJacobian(it->position(),js);
> --
> 1.9.0
>
> From 0a18a5413346a447e054f2ac9c21c525ab6be04b Mon Sep 17 00:00:00 2001
> From: Oliver Sander <sander at igpm.rwth-aachen.de>
> Date: Sat, 29 Mar 2014 11:44:07 +0100
> Subject: [PATCH 5/5] Move computation of quadrature point factor out of inner
> loop
>
> It only depends on the quadrature point and its position. Previously
> it was to deep in a loop structure, and hence computed too many times.
> ---
> dune/pdelab/localoperator/linearelasticity.hh | 12 ++++++------
> 1 file changed, 6 insertions(+), 6 deletions(-)
>
> diff --git a/dune/pdelab/localoperator/linearelasticity.hh b/dune/pdelab/localoperator/linearelasticity.hh
> index 8e001f1..fc9da83 100644
> --- a/dune/pdelab/localoperator/linearelasticity.hh
> +++ b/dune/pdelab/localoperator/linearelasticity.hh
> @@ -103,11 +103,11 @@ namespace Dune {
> RF mu = param_.mu(eg.entity(),it->position());
> RF lambda = param_.lambda(eg.entity(),it->position());
>
> + // geometric weight
> + RF factor = it->weight() * eg.geometry().integrationElement(it->position());
> +
> for(int d=0; d<dim; ++d)
> {
> - // geometric weight
> - RF factor = it->weight() * eg.geometry().integrationElement(it->position());
> -
> for (size_type i=0; i<lfsu.child(0).size(); i++)
> {
> for (int k=0; k<dim; k++)
> @@ -181,6 +181,9 @@ namespace Dune {
> RF mu = param_.mu(eg.entity(),it->position());
> RF lambda = param_.lambda(eg.entity(),it->position());
>
> + // geometric weight
> + RF factor = it->weight() * eg.geometry().integrationElement(it->position());
> +
> for(int d=0; d<dim; ++d)
> {
> const LFSU & lfsu = lfsu_hat.child(d);
> @@ -192,9 +195,6 @@ namespace Dune {
> gradu.axpy(x(lfsu,i),gradphi[i]);
> }
>
> - // geometric weight
> - RF factor = it->weight() * eg.geometry().integrationElement(it->position());
> -
> for (size_type i=0; i<lfsv.child(d).size(); i++)
> {
> for (int k=0; k<dim; k++)
> --
> 1.9.0
>
> _______________________________________________
> 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