functional_aux.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2002-2023 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef FUNCTIONAL_AUX_HH
#define FUNCTIONAL_AUX_HH
 
#include <type_traits>
#include <utility>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include "fem/functionspace.hh"
namespace Kaskade
{
  template <ProblemType Type>
  class FunctionalBase
  {
  public:
 
    static constexpr ProblemType type = Type;
 
    
    // TODO: docme
    template <class T>
    bool inDomain(T const&) const
    {
      return true;
    }
    
    template <int row>
    struct D1
    {
      static constexpr bool present = true;
      
      static constexpr int derivatives = 1;
    };
 
    template <int row, int col>
    struct D2
    {
      static constexpr bool present = type==WeakFormulation || row>=col;
      static constexpr bool symmetric = type==VariationalFunctional && row==col;
      
      static constexpr bool makePositive = false;
 
      static constexpr bool lumped = false;
 
      static constexpr bool constant = false;
 
      static constexpr int derivatives = 1;
    };
 
    template <int row, int col>
    struct B2
    {
      static bool const present = true;
      static bool const symmetric = false;
      static bool const lumped = false;
      static bool const constant = false;
    };
    
    template <class Cell>
    int integrationOrder(Cell const& /* cell */, int shapeFunctionOrder, bool /* boundary */) const 
    {
      return 2*shapeFunctionOrder;
    }
    
    template <class Face>
    bool considerFace(Face const& ) const
    {
      return true;
    }
 
    template <int /*row*/, int /*col*/>
    double makePositiveThreshold() const
    {
      return -10*std::numeric_limits<double>::epsilon();
    }
  };
 
  // ----------------------------------------------------------------------------------------------
 
  template <class TestVars, int row>
  using D1Result = Dune::FieldVector<typename TestVars::Scalar,
                                     TestVars::template Components<row>::m>;
 
 
  template <class TestVars, int row, class AnsatzVars, int col>
  using D2Result = Dune::FieldMatrix<typename TestVars::Scalar,
                                     TestVars::template Components<row>::m,
                                     AnsatzVars::template Components<col>::m>;
 
  // ----------------------------------------------------------------------------------------------
 
  namespace Functional_Aux_Detail
  {
    template <class Functional, int n>
    struct checkSymmetry
    {
      static constexpr bool ok = (!Functional::template D2<n,n>::present)
                                 || Functional::template D2<n,n>::symmetric;
      static constexpr bool pass = Functional::type != VariationalFunctional
                                   || (ok && checkSymmetry<Functional,n-1>::pass);
    };
 
    template <class Functional>
    struct checkSymmetry<Functional,-1>
    {
      static constexpr bool pass = true;
    };
  }
 
  template <class Functional, int n>
  constexpr bool symmetryCheck = Functional_Aux_Detail::checkSymmetry<Functional,n>::pass;
 
 
  // ----------------------------------------------------------------------------------------------
 
 
  namespace Functional_Aux_Detail
  {
    template <int a, int b>
    struct LessThan
    {
      static constexpr bool value = a<b;
    };
 
    template <int a, int b>
    struct Equals
    {
      static constexpr bool value = a==b;
    };
  }
 
  template <class AnsatzVars, class TestVars>
  class TrivialCacheBase
  {
  public:
    using Scalar = typename AnsatzVars::Scalar;
    
    template <class Entity>
    void moveTo(Entity const&) {}
    
    template <class Position, class Evaluators>
    void evaluateAt(Position const&, Evaluators const&) {}
    
    Scalar d0()
    {
      return 0;
    }
    
    template <int row, int dim>
    Dune::FieldVector<Scalar,TestVars::template components<row>> d1(VariationalArg<Scalar,dim> const& arg) const
    {
      return 0;
    }
    
    template <int row, int col, int dim>
    Dune::FieldMatrix <Scalar,TestVars::template components<row>,AnsatzVars::template components<row>>
    d2(VariationalArg<Scalar,dim> const& , VariationalArg<Scalar,dim> const& ) const
    {
      return 0;
    }
    
  };
 
  template <class Functional, class Cache, class ScalarT = typename Functional::Scalar>
  class CacheBase: public TrivialCacheBase<typename Functional::AnsatzVars,typename Functional::TestVars>
  {
  public:
    using Scalar = ScalarT;
    using AnsatzVars = typename Functional::AnsatzVars;
    using TestVars = typename Functional::TestVars;
 
  private:
    template <int row>          using Vector = Dune::FieldVector<Scalar, TestVars::template Components<row>::m>;
    template <int row, int col> using Matrix = Dune::FieldMatrix<Scalar, TestVars::template Components<row>::m, AnsatzVars::template Components<col>::m>;
    template <int row>          using TestComponents =   std::integral_constant<int,TestVars::template Components<row>::m>;
    template <int row>          using AnsatzComponents = std::integral_constant<int,AnsatzVars::template Components<row>::m>;
    
  public:
    template <int row, int dimw>
    Scalar d1_local(VariationalArg<Scalar,dimw,TestComponents<row>::value> const& arg) const
    {
      return static_cast<Cache const*>(this)->template d1_impl<row>(arg);
    }
 
    template <int row, int dimw>
    Dune::FieldVector<Scalar,1> d1(VariationalArg<Scalar,dimw,TestComponents<row>::value> const& arg) const
    {
      return Dune::FieldVector<Scalar,1>(d1_local<row>(arg));
    }  
 
    template <int row, int dimw, class enable = typename std::enable_if_t<Functional_Aux_Detail::LessThan<1,TestComponents<row>::value>::value>>
    Vector<row> d1(VariationalArg<Scalar,dimw> const& scalarArg) const
    {
      Vector<row> result(0);
      VariationalArg<Scalar,dimw,TestComponents<row>::value> arg;
 
      for(int i=0; i<TestComponents<row>::value; ++i)
      {
        arg.value = 0; arg.value[i] = scalarArg.value[0];
        arg.derivative = 0; arg.derivative[i] = scalarArg.derivative[0];
 
        result[i] = d1_local<row>(arg);
      }
 
      return result;
    }
 
    template <int row, int col, int dim>
    Scalar d2_local(VariationalArg<Scalar,dim,TestComponents<row>::value> const& arg1, VariationalArg<Scalar,dim,AnsatzComponents<col>::value> const& arg2) const 
    {
      return static_cast<Cache const*>(this)->template d2_impl<row,col>(arg1,arg2);
    }
 
    template <int row, int col, int dim,
    class enable1 = typename std::enable_if_t<Functional_Aux_Detail::LessThan<1,TestComponents<row>::value>::value>,
    class enable2 = typename std::enable_if_t<Functional_Aux_Detail::LessThan<1,AnsatzComponents<row>::value>::value>
    >
    Dune::FieldMatrix<Scalar,1,1> d2(VariationalArg<Scalar,dim,TestComponents<row>::value> const& arg1, VariationalArg<Scalar,dim,AnsatzComponents<col>::value> const& arg2) const 
    {
      return d2_local<row,col>(arg1,arg2);
    }
 
    template <int row, int col, int dim, int m2,
              class enable1 = typename std::enable_if<Functional_Aux_Detail::LessThan<1,TestComponents<row>::value>::value>::type,
              class enable2 = typename std::enable_if<Functional_Aux_Detail::Equals<m2,AnsatzComponents<col>::value>::value>::type,
              class enable3 = typename std::enable_if<Functional_Aux_Detail::LessThan<1,AnsatzComponents<col>::value>::value>::type
    >
    Vector<row> d2(VariationalArg<Scalar,dim> const& scalarArg, VariationalArg<Scalar,dim,m2> const& arg2) const
    {
      Vector<row> result(0);
      VariationalArg<Scalar,dim,TestComponents<row>::value> arg1;
      for(int i=0; i<AnsatzComponents<row>::value; ++i)
      {
        arg1.value = 0; arg1.value[i] = scalarArg.value[0];
        arg1.derivative = 0; arg1.derivative[i] = scalarArg.derivative[0];
 
        result[i] = d2_local<row,col>(arg1,arg2);
      }
 
      return result;
    }
 
    template <int row, int col, int dim, int m1,
              class enable1 = typename std::enable_if<Functional_Aux_Detail::Equals<m1,TestComponents<row>::value>::value>::type,
              class enable2 = typename std::enable_if<Functional_Aux_Detail::LessThan<1,TestComponents<row>::value>::value>::type,
              class enable3 = typename std::enable_if<Functional_Aux_Detail::LessThan<1,AnsatzComponents<col>::value>::value>::type
             >
    Vector<row> d2(VariationalArg<Scalar,dim,m1> const& arg1, VariationalArg<Scalar,dim> const& scalarArg) const
    {
      Vector<row> result(0);
      VariationalArg<Scalar,dim,AnsatzComponents<row>::value> arg2;
      for(int i=0; i<AnsatzComponents<row>::value; ++i)
      {
        arg2.value = 0; arg2.value[i] = scalarArg.value[0];
        arg2.derivative = 0; arg2.derivative[i] = scalarArg.derivative[0];
 
        result[i] = d2_local<row,col>(arg1,arg2);
      }
 
      return result;
    }
 
    template <int row, int col, int dim>
    Matrix<row,col> d2(VariationalArg<Scalar,dim> const& scalarArg1, VariationalArg<Scalar,dim> const& scalarArg2) const
    {
      Matrix<row,col> result(0);
      VariationalArg<Scalar,dim,TestComponents<row>::value> arg1;
      VariationalArg<Scalar,dim,AnsatzComponents<col>::value> arg2;
      for(int i=0; i<TestComponents<row>::value; ++i)
      {
        arg1.value = 0; arg1.value[i] = scalarArg1.value[0];
        arg1.derivative = 0; arg1.derivative[i] = scalarArg1.derivative[0];
        for(int j=0; j<AnsatzComponents<col>::value; ++j)
        {
          arg2.value = 0; arg2.value[j] = scalarArg2.value[0];
          arg2.derivative = 0; arg2.derivative[j] = scalarArg2.derivative[0];
 
          result[i][j] = d2_local<row,col>(arg1,arg2);
        }
      }
 
      return result;
    }
 
 
    template <int row, int col, int dim>
    Scalar b2_local(VariationalArg<Scalar,dim,TestComponents<row>::value> const& arg1, VariationalArg<Scalar,dim,AnsatzComponents<col>::value> const& arg2) const
    {
      return static_cast<Cache const*>(this)->template b2_impl<row,col>(arg1,arg2);
    }
 
    template <int row, int col, int dim,
    class enable1 = typename std::enable_if<Functional_Aux_Detail::LessThan<1,TestComponents<row>::value>::value>::type,
    class enable2 = typename std::enable_if<Functional_Aux_Detail::LessThan<1,AnsatzComponents<row>::value>::value>::type
    >
    Dune::FieldMatrix<Scalar,1,1> b2(VariationalArg<Scalar,dim,TestComponents<row>::value> const& arg1, VariationalArg<Scalar,dim,AnsatzComponents<col>::value> const& arg2) const
    {
      return b2_local<row,col>(arg1,arg2);
    }
 
    template <int row, int col, int dim, int m2,
    class enable1 = typename std::enable_if<Functional_Aux_Detail::LessThan<1,TestComponents<row>::value>::value>::type,
    class enable2 = typename std::enable_if<Functional_Aux_Detail::Equals<m2,AnsatzComponents<col>::value>::value>::type,
    class enable3 = typename std::enable_if<Functional_Aux_Detail::LessThan<1,AnsatzComponents<col>::value>::value>::type
    >
    Vector<row> b2(VariationalArg<Scalar,dim> const& scalarArg, VariationalArg<Scalar,dim,m2> const& arg2) const
    {
      Vector<row> result(0);
      VariationalArg<Scalar,dim,TestComponents<row>::value> arg1;
      for(int i=0; i<AnsatzComponents<row>::value; ++i)
      {
        arg1.value = 0; arg1.value[i] = scalarArg.value[0];
        arg1.derivative = 0; arg1.derivative[i] = scalarArg.derivative[0];
 
        result[i] = b2_local<row,col>(arg1,arg2);
      }
 
      return result;
    }
 
    template <int row, int col, int dim, int m1,
    class enable1 = typename std::enable_if<Functional_Aux_Detail::Equals<m1,TestComponents<row>::value>::value>::type,
    class enable2 = typename std::enable_if<Functional_Aux_Detail::LessThan<1,TestComponents<row>::value>::value>::type,
    class enable3 = typename std::enable_if<Functional_Aux_Detail::LessThan<1,AnsatzComponents<col>::value>::value>::type
    >
    Vector<row> b2(VariationalArg<Scalar,dim,m1> const& arg1, VariationalArg<Scalar,dim> const& scalarArg) const
    {
      Vector<row> result(0);
      VariationalArg<Scalar,dim,AnsatzComponents<row>::value> arg2;
      for(int i=0; i<AnsatzComponents<row>::value; ++i)
      {
        arg2.value = 0; arg2.value[i] = scalarArg.value[0];
        arg2.derivative = 0; arg2.derivative[i] = scalarArg.derivative[0];
 
        result[i] = b2_local<row,col>(arg1,arg2);
      }
 
      return result;
    }
 
    template <int row, int col, int dim>
    Matrix<row,col> b2(VariationalArg<Scalar,dim> const& scalarArg1, VariationalArg<Scalar,dim> const& scalarArg2) const
    {
      Matrix<row,col> result(0);
      VariationalArg<Scalar,dim,TestComponents<row>::value> arg1;
      VariationalArg<Scalar,dim,AnsatzComponents<col>::value> arg2;
      for(int i=0; i<TestComponents<row>::value; ++i)
      {
        arg1.value = 0; arg1.value[i] = scalarArg1.value[0];
        arg1.derivative = 0; arg1.derivative[i] = scalarArg1.derivative[0];
        for(int j=0; j<AnsatzComponents<col>::value; ++j)
        {
          arg2.value = 0; arg2.value[j] = scalarArg2.value[0];
          arg2.derivative = 0; arg2.derivative[j] = scalarArg2.derivative[0];
 
          result[i][j] = b2_local<row,col>(arg1,arg2);
        }
      }
 
      return result;
    }
  };
  
  // ----------------------------------------------------------------------------
 
  namespace Functional_Aux_Detail
  {
    // helper function that detects if InnerBoundaryCache is declared in Functional
    // If it exists, calling this with a nullptr argument selects this more specific template,
    // otherwise the less specific version below.
    template <class Functional>
    constexpr std::true_type hasInnerBoundaryCacheImpl(typename Functional::InnerBoundaryCache*);
 
    template <class Functional>
    constexpr std::false_type hasInnerBoundaryCacheImpl(...);
  }
 
  template <class Functional>
  constexpr bool hasInnerBoundaryCache =
  decltype(Functional_Aux_Detail::hasInnerBoundaryCacheImpl<Functional>(nullptr))::value;
 
  
 
  // ----------------------------------------------------------------------------
 
  namespace Functional_Aux_Detail
  {
    template <class Functional, class Linearization, bool hasInnerBoundaryCache> 
    struct InnerBoundaryPolicyImpl
    {
    };
 
    template <class Functional, class Linearization>
    struct InnerBoundaryPolicyImpl<Functional,Linearization,true>
    {
      using InnerBoundaryCache = typename Functional::InnerBoundaryCache;
 
      InnerBoundaryCache createInnerBoundaryCache(int flags) const 
      { 
        return InnerBoundaryCache(static_cast<Linearization const*>(this)->getFunctional(),
                                  static_cast<Linearization const*>(this)->getOrigin(),flags);
      }
 
      template <class Face>
      bool considerFace(Face const& face) const
      {
        return static_cast<Linearization const*>(this)->getFunctional().considerFace(face);
      }
    };
 
    template <class Functional,class Linearization> 
    using InnerBoundaryPolicy = InnerBoundaryPolicyImpl<Functional,Linearization,hasInnerBoundaryCache<Functional>>;
  }
 
  template<class Functional_>
  class LinearizationAt : public Functional_Aux_Detail::InnerBoundaryPolicy<Functional_,LinearizationAt<Functional_>>
  {
  public:
    typedef Functional_ Functional;
    static ProblemType const type = Functional::type;
 
    typedef typename Functional::Scalar Scalar;
 
    typedef typename Functional::AnsatzVars    AnsatzVars;
    typedef typename Functional::TestVars      TestVars;
    typedef typename Functional::OriginVars    OriginVars;
 
    typedef typename Functional::DomainCache   DomainCache;
    typedef typename Functional::BoundaryCache BoundaryCache;
 
    typedef typename Functional::OriginVars::VariableSet PointOfLinearization;
 
    LinearizationAt(Functional const & f_, PointOfLinearization const& u_) 
    : f(f_)
    , u(u_)
    { }
    
    DomainCache createDomainCache(int flags) const 
    { 
      return DomainCache(f,u,flags); 
    }
    
    BoundaryCache createBoundaryCache(int flags) const 
    { 
      return BoundaryCache(f,u,flags); 
    }
 
    template <int row>
    struct D1: public Functional::template D1<row> {};
 
    template <int row, int col>
    struct D2: public Functional::template D2<row,col> {};
 
    template <class Cell>
    int integrationOrder(Cell const& cell, int shapeFunctionOrder, bool boundary) const
    {
      return f.integrationOrder(cell,shapeFunctionOrder,boundary);
    }
 
    template <int row, int col>
    auto makePositiveThreshold() const
    {
      return f.template makePositiveThreshold<row,col>();
    }
 
    PointOfLinearization const& getOrigin() const 
    {
      return u;
    }
 
    Functional const& getFunctional() const 
    {
      return f;
    }
    
    int derivatives() const 
    { 
      return f.derivatives(); 
    }
 
    protected:
    Functional const& f;
    PointOfLinearization const& u;
  };
 
 
 
  template<class Functional>
  LinearizationAt<Functional> linearization(Functional const& f, typename Functional::OriginVars::VariableSet const & u)
  {
    return LinearizationAt<Functional>(f,u);
  }
 
  template<class Functional_>
  class LinearizationDifferenceAt : public Functional_Aux_Detail::InnerBoundaryPolicy<Functional_,LinearizationDifferenceAt<Functional_> >
  {
  public:
    typedef Functional_ Functional;
    static ProblemType const type = Functional::type;
 
    typedef typename Functional::Scalar Scalar;
 
    typedef typename Functional::AnsatzVars    AnsatzVars;
    typedef typename Functional::TestVars      TestVars;
    typedef typename Functional::OriginVars    OriginVars;
 
 
    typedef typename Functional::OriginVars::VariableSet PointOfLinearization;
 
    LinearizationDifferenceAt(Functional const & f_,PointOfLinearization const& u1_, PointOfLinearization const& u2_) 
    : f(f_), u1(u1_), u2(u2_) {};
    
    class DomainCache
    {
    public:
      DomainCache(Functional const& f, PointOfLinearization const& u1, PointOfLinearization const& u2, int flags)
      : dc1(f,u1,flags), dc2(f,u2,flags) {}
      
      template <class Cell>
      void moveTo(Cell const& entity) { dc1.moveTo(entity); dc2.moveTo(entity); }
 
      template <class Position, class Evaluators>
      void evaluateAt(Position const& x, Evaluators const& evaluators)  { dc1.evaluateAt(x,evaluators); dc2.evaluateAt(x,evaluators); }
 
      Scalar d0() const { return dc1.d0()-dc2.d0(); }
    
      template<int row, int dim> 
      Dune::FieldVector<Scalar, TestVars::template Components<row>::m>
      d1 (VariationalArg<Scalar,dim> const& arg) const
      {
        return dc1.template d1<row>(arg)-dc2.template d1<row>(arg);
      }
 
      template<int row, int col, int dim> 
      Dune::FieldMatrix<Scalar, TestVars::template Components<row>::m, AnsatzVars::template Components<col>::m>
      d2 (VariationalArg<Scalar,dim> const &arg1, VariationalArg<Scalar,dim> const &arg2) const
      {
        return dc1.template d2<row,col>(arg1,arg2)-dc2.template d2<row,col>(arg1,arg2);
      }
      
      
    private:
      typename Functional::DomainCache dc1, dc2;
    };
    
    class BoundaryCache
    {
    public:
      BoundaryCache(Functional const& f, PointOfLinearization const& u1, PointOfLinearization const& u2, int flags)
      : bc1(f,u1,flags), bc2(f,u2,flags) {}
      
      template <class FaceIterator>
      void moveTo(FaceIterator const& entity) { bc1.moveTo(entity); bc2.moveTo(entity); }
 
      template <class Position, class Evaluators>
      void evaluateAt(Position const& x, Evaluators const& evaluators) { bc1.evaluateAt(x,evaluators); bc2.evaluateAt(x,evaluators); }
 
      Scalar d0() const { return bc1.d0()-bc2.d0(); }
    
      template<int row, int dim> 
      Dune::FieldVector<Scalar, TestVars::template Components<row>::m> 
      d1 (VariationalArg<Scalar,dim> const& arg) const
      {
        return bc1.template d1<row>(arg)-bc2.template d1<row>(arg);
      }
 
      template<int row, int col, int dim> 
      Dune::FieldMatrix<Scalar, TestVars::template Components<row>::m, AnsatzVars::template Components<col>::m>
      d2 (VariationalArg<Scalar,dim> const &arg1, VariationalArg<Scalar,dim> const &arg2) const
      {
        return bc1.template d2<row,col>(arg1,arg2)-bc2.template d2<row,col>(arg1,arg2);
      }
      
    private:
      typename Functional::BoundaryCache bc1, bc2;
    };
 
    DomainCache createDomainCache(int flags) const 
    {
      return DomainCache(f,u1,u2,flags);
    }
    
    BoundaryCache createBoundaryCache(int flags) const 
    {
      return BoundaryCache(f,u1,u2,flags);
    }
 
    template <int row>
    struct D1: public Functional::template D1<row> {};
 
    template <int row, int col>
    struct D2: public Functional::template D2<row,col> {};
 
    template <class Cell>
    int integrationOrder(Cell const& cell, int shapeFunctionOrder, bool boundary) const
    {
      return f.template integrationOrder<Cell>(cell, shapeFunctionOrder, boundary);
    }
 
    Functional const& getFunctional() const {return f;}
 
  protected:
    Functional const& f;
    PointOfLinearization const& u1;
    PointOfLinearization const& u2;
  };
 
 
 
  template<class Functional>
  LinearizationDifferenceAt<Functional> linearizationDifference(Functional const& f, 
                                                                typename Functional::OriginVars::VariableSet const & u1, 
                                                                typename Functional::OriginVars::VariableSet const & u2)
  {
    return LinearizationDifferenceAt<Functional>(f,u1,u2);
  }
 
  //---------------------------------------------------------------------
 
  template <class Functional>
  class SemiLinearizationAt: public LinearizationAt<Functional>
  {
  public:
    SemiLinearizationAt(Functional const& f_,
                        typename Functional::AnsatzVars::VariableSet const& u_,
                        typename Functional::AnsatzVars::VariableSet const& uJ_,
                        typename Functional::AnsatzVars::VariableSet const& du_):
          LinearizationAt<Functional>(f_,u_), uJ(uJ_), du(du_) {}
 
    typename Functional::DomainCache createDomainCache(int flags) const 
    {
      return typename Functional::DomainCache(this->f,this->u,uJ,du,flags);
    }
 
    typename Functional::BoundaryCache createBoundaryCache(int flags) const 
    {
      return typename Functional::BoundaryCache(this->f,this->u,uJ,du,flags);
    }
    
  private:
    typename Functional::AnsatzVars::VariableSet const& uJ;
    typename Functional::AnsatzVars::VariableSet const& du;
  };
 
  template <class Functional>
  class SemiLinearizationAtInner: public LinearizationAt<Functional>
  {
  public:
    SemiLinearizationAtInner(Functional const& f_,
                        typename Functional::AnsatzVars::VariableSet const& u_,
                        typename Functional::AnsatzVars::VariableSet const& uJ_,
                        typename Functional::AnsatzVars::VariableSet const& du_):
          LinearizationAt<Functional>(f_,u_), uJ(uJ_), du(du_) {}
 
    typename Functional::DomainCache createDomainCache(int flags) const 
    {
      return typename Functional::DomainCache(this->f,this->u,uJ,du,flags);
    }
 
    typename Functional::BoundaryCache createBoundaryCache(int flags) const 
    {
      return typename Functional::BoundaryCache(this->f,this->u,uJ,du,flags);
    }
    
    typename Functional::InnerBoundaryCache createInnerBoundaryCache(int flags) const 
    {
      return typename Functional::InnerBoundaryCache(this->f,this->u,uJ,du,flags);
    }
    
  private:
    typename Functional::AnsatzVars::VariableSet const& uJ;
    typename Functional::AnsatzVars::VariableSet const& du;
  };
 
  template<class Functional>
  SemiLinearizationAt<Functional> semiLinearization(Functional const& f,
      typename Functional::AnsatzVars::VariableSet const& u,
      typename Functional::AnsatzVars::VariableSet const& uJ,
      typename Functional::AnsatzVars::VariableSet const& du)
  {
    return SemiLinearizationAt<Functional>(f,u,uJ,du);
  }
 
  //---------------------------------------------------------------------------------------------
  //---------------------------------------------------------------------------------------------
  
  
} /* end of namespace Kaskade */
 
#endif