[dune-pdelab] Patches for local operators

Oliver Sander sander at igpm.rwth-aachen.de
Sat Apr 5 16:18:59 CEST 2014


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
> 

-------------- next part --------------
A non-text attachment was scrubbed...
Name: 0001-Use-the-new-Laplace-assembler-in-laplace.hh-to-imple.patch
Type: text/x-diff
Size: 3645 bytes
Desc: not available
URL: <https://lists.dune-project.org/pipermail/dune-pdelab/attachments/20140405/d528aa06/attachment.patch>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 0002-Implement-the-matrix-assemble-i.e.-jacobian_volume-i.patch
Type: text/x-diff
Size: 2079 bytes
Desc: not available
URL: <https://lists.dune-project.org/pipermail/dune-pdelab/attachments/20140405/d528aa06/attachment-0001.patch>
-------------- 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/d528aa06/attachment.sig>


More information about the dune-pdelab mailing list