variables.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-2024 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef MULTIVARIATIONALUTIL_HH
#define MULTIVARIATIONALUTIL_HH
 
#include <algorithm>
#include <array>
#include <cassert>
#include <iostream>
 
#include <boost/fusion/algorithm.hpp>
#include <boost/fusion/include/fold.hpp>
#include <boost/fusion/include/make_vector.hpp>
#include <boost/fusion/include/size.hpp>
#include <boost/fusion/sequence.hpp>
#include <boost/mpl/int.hpp>
#include <boost/type_traits/remove_reference.hpp>
#include <boost/type_traits/remove_const.hpp>
 
 
#include "fem/fixdune.hh"
#include "fem/functionspace.hh"
#include "fem/linearspace.hh"
#include "io/iobase.hh" // for paddedString
#include "linalg/crsutil.hh"
 
//---------------------------------------------------------------------
 
namespace Kaskade
{
  template <int spaceId, int components, int Id=-1>
  struct VariableDescription
  {
    static int const id = Id;
    static int const m = components;
    static int const spaceIndex = spaceId;
  };
 
 
  template <int spaceIndex>
  struct SpaceIndex : public VariableDescription<spaceIndex,-1,-1> {};
 
  template <int components>
  struct Components : public VariableDescription<-1,components,-1> {};
 
  template <int id>
  struct [[deprecated]] VariableId : public VariableDescription<-1,-1,id> {};
 
 
  template <typename A, typename B, typename C = VariableDescription<-1, -1, -1>>
  struct Variable 
  {
    static int const spaceIndex = A::spaceIndex>=0? A::spaceIndex : B::spaceIndex>=0? B::spaceIndex : C::spaceIndex;
    static int const m  = A::m>=0? A::m : B::m>=0? B::m : C::m;
 
    static_assert(spaceIndex>=0,"Space index has to be nonnegative (space indices cover a contiguous range from 0 to k-1).");
    static_assert(m>=0,"Number of components has to be nonnegative (in fact, 0 makes rarely sense as well).");
 
    Variable(std::string const& name_)
    : name(name_) {}
    
    Variable() = default;
    
    std::string name;
  };
  
 
  //---------------------------------------------------------------------
 
  template <class Spaces>
  struct ConstructElement
  {
    ConstructElement(Spaces const& spaces_):
      spaces(spaces_)
    {}
 
    template <typename Arg> struct result {};
 
    template <class Variable>
    struct result<ConstructElement<Spaces>(Variable)>
    {
      typedef typename boost::remove_reference<Variable>::type Var;
      typedef typename SpaceType<Spaces,Var::spaceIndex>::type Space;
      typedef typename Space::template Element<Var::m>::type type;
    };
 
    template <class Variable>
    typename result<ConstructElement<Spaces>(Variable)>::type operator()(Variable const& variable) const
    {
      return typename result<ConstructElement<Spaces>(Variable)>::type(*boost::fusion::at_c<Variable::spaceIndex>(spaces));
    }
 
  private:
    Spaces const& spaces;
  };
 
  template <class Spaces>
  struct ConstructCoefficientVector
  {
    ConstructCoefficientVector(Spaces const& spaces_):
      spaces(spaces_)
    {}
 
    template <typename Arg> struct result {};
 
    template <class Variable>
    struct result<ConstructCoefficientVector<Spaces>(Variable)>
    {
      typedef typename boost::remove_reference<Variable>::type Var;
      typedef typename SpaceType<Spaces,Var::spaceIndex>::type Space;
      typedef typename Space::template Element<Var::m>::type::StorageType StorageType;
 
      typedef Dune::BlockVector<typename StorageType::block_type> type;
    };
 
    template <class Variable>
    typename result<ConstructCoefficientVector<Spaces>(Variable)>::type operator()(Variable const& variable) const
    {
      return typename result<ConstructCoefficientVector<Spaces>(Variable)>::type(boost::fusion::at_c<Variable::spaceIndex>(spaces)->degreesOfFreedom());
    }
 
  private:
    Spaces const& spaces;
  };
 
  //---------------------------------------------------------------------
  
  namespace Variables_Detail 
  {
    using namespace boost::fusion;
    
    // A functor constructing representations of consecutive subsets of variables
    template <class Variables, class RepresentationConstructor,
              int first=0, int last=result_of::size<Variables>::type::value>
    class VariableRangeCreator
    {
      typedef typename result_of::begin<Variables>::type       Begin;
      typedef typename result_of::advance_c<Begin,first>::type First;
      typedef typename result_of::advance_c<Begin,last>::type  Last;
      
    public:
      typedef typename boost::fusion::transform_view<iterator_range<First,Last> const,
                                                     RepresentationConstructor> View;
      using type = typename result_of::as_vector<View>::type;
      
      static View apply(RepresentationConstructor const& c)
      {
        using namespace boost::fusion;
        Variables vars;
        return transform(iterator_range<First,Last>(advance_c<first>(begin(vars)),advance_c<last>(begin(vars))),c);
      }
    };
 
    // Metafunction class for boost mpl transform. It takes a variable description and an id,
    // and overwrites the variable id by the given id.
    struct AddIdToVariableDescription 
    {
      template <class Vd, class Int>
      struct apply 
      {
        using type = Kaskade::VariableDescription<Vd::spaceIndex, Vd::m, Int::value>;
      };
    };
 
    // Given a type of a sequence of VariableDescriptions a new sequence type of VariableDescriptions is created 
    // where the id of each VariableDescription is the position in the sequence
    template <typename IncompleteVariableDescriptions>
    struct ImplicitIdVariableDescriptions 
    {
      using size = typename result_of::size<IncompleteVariableDescriptions>::type;
      using type = typename result_of::as_vector<typename boost::mpl::transform<IncompleteVariableDescriptions, 
                                                                                boost::mpl::range_c<int,0,size::value>,
                                                                                AddIdToVariableDescription>::type
                                                 >::type;
    };
    
    // Given a sequence of variables and a sequence of spaces, returns the minimal
    // continuity of all variables.
    template <typename Variables, typename Spaces>
    class Continuity
    {
      struct Min
      {
        template <typename C, typename Var>
        auto operator()(C,Var) const 
        {
          using SpacePtr = typename result_of::value_at_c<Spaces,Var::spaceIndex>::type;
          constexpr int m = std::remove_pointer_t<SpacePtr>::continuity;
          return typename boost::mpl::min<C,boost::mpl::int_<m>>::type();
        }
      };
      
    public:
      static int const value = result_of::fold<Variables,boost::mpl::int_<1000>,Min>::type::value;
    };
    
  }
  //---------------------------------------------------------------------
 
  template <class VSDescriptions>
  class VariableSet: public LinearProductSpace<typename VSDescriptions::Scalar, 
                                               typename VSDescriptions::RepresentationData>
  {
  public:
    typedef VSDescriptions Descriptions;
 
    typedef typename VSDescriptions::RepresentationData Functions;
 
  private:
    typedef VariableSet<Descriptions> Self;
    typedef typename SpaceType<typename Descriptions::Spaces,0>::type::Scalar Scalar;
    typedef LinearProductSpace<Scalar,Functions> Base;
 
    typedef typename Descriptions::Grid Grid;
    using GridView = typename Descriptions::GridView;
 
  public:
    typedef LinearProductSpace<Scalar,Functions> LinearSpace;
 
    VariableSet(Self const& vs):
      Base(vs.data), descriptions(vs.descriptions)
    {}
 
    explicit VariableSet(Descriptions const& d):
        Base(Variables_Detail::VariableRangeCreator<typename Descriptions::Variables,ConstructElement<typename Descriptions::Spaces> >::apply(
          ConstructElement<typename Descriptions::Spaces>(d.spaces))), 
        descriptions(d)
    {}
 
    // TODO: The reference to the variable set description is not modified. Danger is,
    // that the assigned-to function space elements are now elements from different spaces
    // than before (and stated by the variable set description). Maybe switch to a (private)
    // pointer member and provide a constant getter?
    Self& operator=(Self const& v)
    {
      if (this!=&v)
        this->data = v.data;
      return *this;
    }
 
    template <class S>
    Self& operator=(S const& s)
    {
      static_cast<Base&>(*this) = s;
      return *this;
    }
    
    template <int i, class Evaluators>
    auto value(Evaluators const& evals) const 
    {
      static_assert(boost::fusion::result_of::size<Evaluators>::type::value==boost::fusion::result_of::size<typename Descriptions::Spaces>::type::value,
                    "Evaluators have to match the space list");
      constexpr int spaceIndex = Descriptions::template spaceIndex<i>;
      return component<i>(*this).value(boost::fusion::at_c<spaceIndex>(evals));
    }
 
    template <int i>
    auto value(Cell<GridView> const& cell, LocalPosition<GridView> const& xi) const
    {
      return component<i>(*this).value(cell,xi);
    }
    
    template <int i>
    auto value(GlobalPosition<GridView> const& x) const
    {
      return component<i>(*this).value(x);
    }
 
    template <int i, class Evaluators>
    auto derivative(Evaluators const& evals) const
    {
      int const spaceIndex = Descriptions::template spaceIndex<i>;
      return component<i>(*this).derivative(boost::fusion::at_c<spaceIndex>(evals));
    }
    
    template <int i>
    auto derivative(Cell<GridView> const& cell, LocalPosition<GridView> const& xi) const
    {
      return component<i>(*this).derivative(cell,xi);
    }
    
    template <int i>
    auto derivative(GlobalPosition<GridView> const& x) const
    {
      return component<i>(*this).derivative(x);
    }
 
    template <int i, class Evaluators>
    auto hessian(Evaluators const& evals) const
    {
      int const spaceIndex = Descriptions::template spaceIndex<i>;
      return component<i>(*this).hessian(boost::fusion::at_c<spaceIndex>(evals));
    }
    
    template <int i>
    auto hessian(Cell<GridView> const& cell, LocalPosition<GridView> const& xi) const
    {
      return component<i>(*this).hessian(cell,xi);
    }
    
    template <int i>
    auto hessian(GlobalPosition<GridView> const& x) const
    {
      return component<i>(*this).hessian(x);
    }
    
 
    Descriptions const& descriptions;
  };
 
  //---------------------------------------------------------------------
  //---------------------------------------------------------------------
 
 
  template <class SpaceList, class VariableList>
  class VariableSetDescription
  {
  public:
    typedef SpaceList Spaces;
    
    using Variables = typename Variables_Detail::ImplicitIdVariableDescriptions<VariableList>::type;
 
    typedef typename Variables_Detail::VariableRangeCreator<Variables,ConstructElement<Spaces>>::type RepresentationData;
 
    typedef Kaskade::VariableSet<VariableSetDescription<Spaces,VariableList>> VariableSet;
    
    typedef typename SpaceType<Spaces,0>::type::Scalar   Scalar;
    typedef typename SpaceType<Spaces,0>::type::Grid     Grid;
    typedef typename SpaceType<Spaces,0>::type::GridView GridView;
    typedef typename SpaceType<Spaces,0>::type::IndexSet IndexSet;
 
    static int const noOfVariables = boost::fusion::result_of::size<Variables>::type::value;
    
    static int const continuity = Variables_Detail::Continuity<Variables,Spaces>::value;
    
    typedef typename boost::fusion::result_of::as_vector<
          typename boost::fusion::result_of::transform<Spaces, GetEvaluators<ShapeFunctionCache<Grid,Scalar>> >::type
          >::type Evaluators;
 
    template <class RAIter>
    VariableSetDescription(Spaces const& spaces_, RAIter nameIt):
      VariableSetDescription(spaces_)
    {
      std::copy(nameIt,nameIt+noOfVariables,names.begin());
    }
 
    VariableSetDescription(Spaces const& spaces_, std::vector<std::string> names_)
    : VariableSetDescription(spaces_)
    {
      assert(names.size()==names_.size());
      std::copy(names_.begin(),names_.end(),names.begin());
    }
 
    template <class RAIter>
    VariableSetDescription(Spaces const& spaces_, RAIter nameIt, RAIter roleIt)
    : VariableSetDescription(spaces_,nameIt)
    {
      std::copy(roleIt,roleIt+noOfVariables,roles.begin());
    }
 
    VariableSetDescription(Spaces const& spaces_, std::vector<std::string> names_, std::vector<std::string> roles_)
    : VariableSetDescription(spaces_,names_)
    {
      std::copy(roles_.begin(),roles_.end(),roles.begin());
    }
 
    VariableSetDescription(Spaces const& spaces_)
    : spaces(spaces_)
    , gridView(boost::fusion::at_c<0>(spaces_)->gridView())
    , indexSet(boost::fusion::at_c<0>(spaces_)->indexSet())
    {}
 
    size_t degreesOfFreedom(int first=0, int last=noOfVariables) const
    {
      return degreesOfFreedom(spaces,first,last);
    }
 
    static size_t degreesOfFreedom(Spaces const& spaces, int first, int last)
    {
      auto dimensions = variableDimensions(spaces);
      assert(0<=first && first<=last && last<=dimensions.size());
      return std::accumulate(begin(dimensions)+first,begin(dimensions)+last,0);
    }      
    
    static std::array<size_t,noOfVariables> variableDimensions(Spaces const& spaces, int first=0, int last=noOfVariables)
    {
      std::array<size_t,noOfVariables> dimensions;
      boost::fusion::for_each(Variables(),ComputeDimension<noOfVariables>(dimensions,spaces));
      return dimensions;
    }
 
    template <int i>
    constexpr static int components = boost::fusion::result_of::value_at_c<Variables,i>::type::m;
 
    template <int idx>
    struct [[deprecated]] Components { static int const m = boost::fusion::result_of::value_at_c<Variables,idx>::type::m; };
 
    template <int idx>
    constexpr static int spaceIndex = boost::fusion::result_of::value_at_c<Variables,idx>::type::spaceIndex;
 
    template <int idx>
    struct [[deprecated]] SpaceIndex { static int const spaceIndex = boost::fusion::result_of::value_at_c<Variables,idx>::type::spaceIndex; };
 
    
    template <int idx>
    using Space = std::remove_pointer_t<typename boost::fusion::result_of::value_at_c<Spaces,SpaceIndex<idx>::spaceIndex>::type>;
    
    template <int first=0, int last=noOfVariables>
    class [[deprecated]] CoefficientVectorRepresentation
    {
      using Creator = Variables_Detail::VariableRangeCreator<Variables,ConstructCoefficientVector<Spaces>,first,last>;
 
    public:
      using type = LinearProductSpace<Scalar,typename Creator::type>;
      
      static typename Creator::View init(Spaces const& spaces) 
      {
        return Creator::apply(ConstructCoefficientVector<Spaces>(spaces));
      }
    };
 
    template <int first=0, int last=noOfVariables>
    using CoefficientVector = typename CoefficientVectorRepresentation<first,last>::type;
    
    template <int first=0, int last=noOfVariables>
    auto zeroCoefficientVector() const 
    {
      using CVec = CoefficientVectorRepresentation<first,last>;
      return typename CVec::type(CVec::init(spaces));
    }
    
    VariableSet variableSet() const 
    {
      return VariableSet(*this);
    }
 
    Spaces             spaces;
    
    GridView const&    gridView;
    
    GridManagerBase<Grid> const& gridManager() const { return boost::fusion::at_c<0>(spaces)->gridManager(); }
 
    IndexSet const&    indexSet;
    
    std::array<std::string,noOfVariables>  names;
 
    std::array<std::string,noOfVariables>  roles;
 
  private:
 
    struct CheckHasEqualGridAndIndexSet
    {
      CheckHasEqualGridAndIndexSet(GridManagerBase<Grid> const* gp_, IndexSet const* isp_):
        gp(gp_), isp(isp_)
      {}
 
      template <class Space>
      void operator()(Space const* space) const {
        assert(&(space->gridManager())== gp);
        assert(&(space->indexSet()) == isp);
      }
 
      GridManagerBase<Grid> const* gp;
      IndexSet const* isp;
    };
 
    template <int nvars>
    struct ComputeDimension
    {
      ComputeDimension(std::array<size_t,nvars>& dims, Spaces const& spaces_):
        dims_(dims), spaces(spaces_)
      {}
 
      template <class Variable>
      void operator()(Variable const& /* v */) const
      {
        assert(dims_.size()>Variable::id);
        dims_[Variable::id] = Variable::m * boost::fusion::at_c<Variable::spaceIndex>(spaces)->degreesOfFreedom();
      }
 
      std::array<size_t,nvars>& dims_;
      Spaces const&             spaces;
    };
 
  };
 
  template <class VarSetDesc, int varIdx>
  constexpr int spaceIndex = VarSetDesc::template spaceIndex<varIdx>;
  
  //---------------------------------------------------------------------
  
  
  template <class ...Spaces>
  auto makeSpaceList(Spaces*... spaces)
  {
    return boost::fusion::make_vector(const_cast<Spaces const*>(spaces)...);
  }
  
  template <class SpaceList, class VariableList>
  VariableSetDescription<SpaceList,VariableList> makeVariableSetDescription(SpaceList const& spaces,
                                                                            VariableList const& vars)
  {
    using namespace boost::fusion;
    
    // Extract the names of the variables.
    std::vector<std::string> names;
    for_each(vars,[&](auto var) { names.push_back(var.name); });
    // If (some) variable names are not specified, assign default ones.
    for (int i=0; i<names.size(); ++i)
      if (names[i].empty())
        names[i] = "var" + paddedString(i,2);
      
    return VariableSetDescription<SpaceList,VariableList>(spaces,names);
  }
  
 
  //---------------------------------------------------------------------
  //---------------------------------------------------------------------
  
 
  template <class Variables, class Functions, class Evaluators, class Method>
  auto evaluateFunctions(Functions const& fs,
                         Evaluators const& eval,
                         Method const& method = Method())
  {
    using namespace boost::fusion;
    
    // Define the functor that evaluates the individual variables.
    auto doEval = [&] (auto const& pair)
    {
      // Extract the space index from the second component of the pair provided by zip.
      static int const sid = std::remove_reference_t<typename result_of::value_at_c<std::remove_reference_t<decltype(pair)>,1>::type>::spaceIndex;
      // Evaluate with evaluator corresponding to the function's space.
      return method(at_c<0>(pair),at_c<sid>(eval));
    };
    
    // We need this as vector type to default construct an object for transform. This way we can submit
    // views (such as joint_view) as template parameter.
    using VarDesc = typename result_of::as_vector<Variables>::type;
    
    // For evaluation, the corresponding evaluator must be used based on the space index of the variables. To provide this,
    // we zip the functions with their variable descriptions to pairs from which we can extract both function and space index.
    return as_vector(transform(zip(fs,VarDesc()),doEval));
  }
  
  template <class VariableSet, class Method>
  auto evaluateVariables(VariableSet const& vars,
                         typename VariableSet::Descriptions::Evaluators const& eval,
                         Method const& method = Method())
  {
    return evaluateFunctions<typename VariableSet::Descriptions::Variables>(vars.data,eval,method);
  }
  
  template <class VariableSet, class Method>
  using EvaluateVariables = decltype(evaluateVariables(std::declval<VariableSet>(),
                                                       std::declval<typename VariableSet::Descriptions::Evaluators>(),
                                                       std::declval<Method>()));
 
 
  constexpr auto valueMethod = [](auto const& f, auto const& eval) { return f.value(eval); };
  
  constexpr auto derivativeMethod = [](auto const& f, auto const& eval) { return f.derivative(eval); };
  
  using ValueMethod = decltype(valueMethod);
  
  using DerivativeMethod = decltype(derivativeMethod);
 
} // end of namespace Kaskade
 
 
 
 
//---------------------------------------------------------------------
//---------------------------------------------------------------------
//---------------------------------------------------------------------
 
#endif