[Dune-devel] [Dune-Commit] [Commit] dune-localfunctions - 13481f3: Implement second derivatives for PK2d shape functions

[Carsten Gräser]

Hi,
I'd like to merge this feature to 2.4, but it seems that this
is incomplete, because we guarantee up to 2nd order partial
derivatives while throwing if 1st order is requested.

Is someone working on the implementation of 1st order partial
derivatives? Otherwise I'll drop this for 2.4.

Best,
Carsten



Am 29.04.2015 um 17:13 schrieb Oliver Sander:
> New commit, appeared at Wed Apr 29 17:13:30 2015 +0200
> as part of the following ref changes:
> 
>     branch refs/heads/master    updated from d127aad -> 13481f3
> 
> Browsable version: http://cgit.dune-project.org/repositories/dune-localfunctions/commit/?id=13481f3bb3d4ab1e033967bd9c0a8f28f3b3fe9a
> 
> ======================================================================
> 
> commit 13481f3bb3d4ab1e033967bd9c0a8f28f3b3fe9a
> Author:     Elisa Friebel <friebel at igpm.rwth-aachen.de>
> AuthorDate: Wed Apr 29 17:10:36 2015 +0200
> Commit:     Oliver Sander <sander at igpm.rwth-aachen.de>
> CommitDate: Wed Apr 29 17:10:36 2015 +0200
> 
>     Implement second derivatives for PK2d shape functions
> 
>  COPYING                                            |  1 +
>  .../localfunctions/lagrange/pk2d/pk2dlocalbasis.hh | 79 +++++++++++++++++++++-
>  2 files changed, 79 insertions(+), 1 deletion(-)
> 
> 
> 
> diff --git a/COPYING b/COPYING
> index 4ec1a66..308b0a1 100644
> --- a/COPYING
> +++ b/COPYING
> @@ -10,6 +10,7 @@ Copyright holders:
>  2008--2015    Christian Engwer
>  2008--2012    Jorrit Fahlke
>  2011--2013    Bernd Flemisch
> +2015          Elisa Friebel
>  2012          Christoph Gersbacher
>  2009--2014    Carsten Gräser
>  2015          Felix Gruber
> diff --git a/dune/localfunctions/lagrange/pk2d/pk2dlocalbasis.hh b/dune/localfunctions/lagrange/pk2d/pk2dlocalbasis.hh
> index 8286e19..0ccea9c 100644
> --- a/dune/localfunctions/lagrange/pk2d/pk2dlocalbasis.hh
> +++ b/dune/localfunctions/lagrange/pk2d/pk2dlocalbasis.hh
> @@ -35,7 +35,7 @@ namespace Dune
>      enum {O = k};
>  
>      typedef LocalBasisTraits<D,2,Dune::FieldVector<D,2>,R,1,Dune::FieldVector<R,1>,
> -        Dune::FieldMatrix<R,1,2> > Traits;
> +        Dune::FieldMatrix<R,1,2>, 2 > Traits;
>  
>      //! \brief Standard constructor
>      Pk2DLocalBasis ()
> @@ -158,6 +158,59 @@ namespace Dune
>  
>      }
>  
> +    //! \brief Evaluate higher derivatives of all shape functions
> +    template<unsigned int order> //order of derivative
> +    inline void evaluate(const std::array<int,order>& directions, //direction of derivative
> +                         const typename Traits::DomainType& in,  //position
> +                         std::vector<typename Traits::RangeType>& out) const //return value
> +    {
> +      out.resize(N);
> +
> +      if (order > Traits::diffOrder)
> +        DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
> +
> +      if (order==0)
> +        evaluateFunction(in, out);
> +      else if (order==1)
> +        DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
> +      else if (order==2)
> +      {
> +        // specialization for k<2, not clear whether that is needed
> +        if (k<2) {
> +          std::fill(out.begin(), out.end(), 0.0);
> +        return;
> +      }
> +
> +      //f = prod_{i} f_i -> dxa dxb f = sum_{i} {dxa f_i sum_{k \neq i} dxb f_k prod_{l \neq k,i} f_l
> +      int n=0;
> +      for (unsigned int j=0; j<=k; j++)
> +        for (unsigned int i=0; i<=k-j; i++, n++)
> +        {
> +          R res = 0.0;
> +
> +          for (unsigned int no1=0; no1 < k; no1++)
> +          {
> +            R factor1 = lagrangianFactorDerivative(directions[0], no1, i, j, in);
> +            for (unsigned int no2=0; no2 < k; no2++)
> +            {
> +              if (no1 == no2)
> +                continue;
> +              R factor2 = factor1*lagrangianFactorDerivative(directions[1], no2, i, j, in);
> +              for (unsigned int no3=0; no3 < k; no3++)
> +              {
> +                if (no3 == no1 || no3 == no2)
> +                  continue;
> +                factor2 *= lagrangianFactor(no3, i, j, in);
> +              }
> +              res += factor2;
> +            }
> +
> +          }
> +          out[n] = res;
> +        }
> +      }
> +    }
> +
>      //! \brief Polynomial order of the shape functions
>      unsigned int order () const
>      {
> @@ -165,6 +218,30 @@ namespace Dune
>      }
>  
>    private:
> +  /** \brief Returns a single Lagrangian factor of l_ij evaluated at x */
> +  typename Traits::RangeType lagrangianFactor(const int no, const int i, const int j, const typename Traits::DomainType& x) const
> +  {
> +    if ( no < i)
> +      return (x[0]-pos[no])/(pos[i]-pos[no]);
> +    if (no < i+j)
> +      return (x[1]-pos[no-i])/(pos[j]-pos[no-i]);
> +    return (pos[no+1]-x[0]-x[1])/(pos[no+1]-pos[i]-pos[j]);
> +  }
> +
> +  /** \brief Returns the derivative of a single Lagrangian factor of l_ij evaluated at x
> +   * \param direction Derive in x-direction if this is 0, otherwise derive in y direction
> +   */
> +  typename Traits::RangeType lagrangianFactorDerivative(const int direction, const int no, const int i, const int j, const typename Traits::DomainType& x) const
> +  {
> +    if ( no < i)
> +      return (direction == 0) ? 1.0/(pos[i]-pos[no]) : 0;
> +
> +    if (no < i+j)
> +      return (direction == 0) ? 0: 1.0/(pos[j]-pos[no-i]);
> +
> +    return -1.0/(pos[no+1]-pos[i]-pos[j]);
> +  }
> +
>      R pos[k+1]; // positions on the interval
>    };

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 473 bytes
Desc: OpenPGP digital signature
URL: <https://lists.dune-project.org/pipermail/dune-devel/attachments/20150505/60565da7/attachment.sig>