[dune-pdelab] Patches for local operators
Oliver Sander
sander at igpm.rwth-aachen.de
Sat Apr 5 14:10:54 CEST 2014
Hi Christian,
thanks for committing. The Laplace class actually only consists of a single
integer, so the constructor doesn't cost you anything. Do you still insist
on making it a private member of Poisson?
Please also commit patches 4 and 5. They contain cleanup/speedup changes
exclusively, and do not change API or ABI.
Thanks,
Oliver
> 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
>
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 534 bytes
Desc: OpenPGP digital signature
URL: <https://lists.dune-project.org/pipermail/dune-pdelab/attachments/20140405/c27e9ff1/attachment.sig>
More information about the dune-pdelab
mailing list