[dune-pdelab] Patches for local operators
Christian Engwer
christian.engwer at uni-muenster.de
Sat Apr 5 18:21:41 CEST 2014
Hi Oli,
Thanks, I'll check and commit when I'm back at my computer.
Christian
On 5. April 2014 16:18:59 MESZ, Oliver Sander <sander at igpm.rwth-aachen.de> wrote:
>Hi Christian,
>here are the two reworked patches to poisson.hh. The Laplace assembler
>is now a data
>member of the Poisson class.
>Cheers,
>Oliver
>
>Am 05.04.2014 14:10, schrieb Oliver Sander:
>> 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
>>>
>>>
>>
>>
>>
>>
>> _______________________________________________
>> dune-pdelab mailing list
>> dune-pdelab at dune-project.org
>> http://lists.dune-project.org/mailman/listinfo/dune-pdelab
>>
>
>
>
>------------------------------------------------------------------------
>
>_______________________________________________
>dune-pdelab mailing list
>dune-pdelab at dune-project.org
>http://lists.dune-project.org/mailman/listinfo/dune-pdelab
More information about the dune-pdelab
mailing list