gridoperator.hh

#ifndef DUNE_PDELAB_GRIDOPERATOR_GRIDOPERATOR_HH
#define DUNE_PDELAB_GRIDOPERATOR_GRIDOPERATOR_HH
 
#include <tuple>
 
#include <dune/common/hybridutilities.hh>
 
#include <dune/pdelab/gridfunctionspace/interpolate.hh>
#include <dune/pdelab/gridoperator/common/borderdofexchanger.hh>
#include <dune/pdelab/gridoperator/common/gridoperatorutilities.hh>
#include <dune/pdelab/gridoperator/default/assembler.hh>
#include <dune/pdelab/gridoperator/default/localassembler.hh>
 
namespace Dune{
  namespace PDELab{
 
    template<typename GFSU, typename GFSV, typename LOP,
             typename MB, typename DF, typename RF, typename JF,
             typename CU=Dune::PDELab::EmptyTransformation,
             typename CV=Dune::PDELab::EmptyTransformation
             >
    class GridOperator
    {
    public:
 
      typedef DefaultAssembler<GFSU,GFSV,CU,CV> Assembler;
 
      using Domain = Dune::PDELab::Backend::Vector<GFSU,DF>;
      using Range = Dune::PDELab::Backend::Vector<GFSV,RF>;
      using Jacobian = Dune::PDELab::Backend::Matrix<MB,Domain,Range,JF>;
 
      typedef typename MB::template Pattern<Jacobian,GFSV,GFSU> Pattern;
 
      typedef DefaultLocalAssembler<
        GridOperator,
        LOP,
        GFSU::Traits::EntitySet::Partitions::partitionIterator() == InteriorBorder_Partition
        > LocalAssembler;
 
      // Fix this as soon as the default Partitions are constexpr
      typedef typename std::conditional<
        GFSU::Traits::EntitySet::Partitions::partitionIterator() == InteriorBorder_Partition,
        NonOverlappingBorderDOFExchanger<GridOperator>,
        OverlappingBorderDOFExchanger<GridOperator>
        >::type BorderDOFExchanger;
 
      typedef Dune::PDELab::GridOperatorTraits
      <GFSU,GFSV,MB,DF,RF,JF,CU,CV,Assembler,LocalAssembler> Traits;
 
      template <typename MFT>
      struct MatrixContainer{
        typedef typename Traits::Jacobian Type;
      };
 
      GridOperator(const GFSU & gfsu_, const CU & cu_, const GFSV & gfsv_, const CV & cv_, LOP & lop_, const MB& mb_ = MB())
        : global_assembler(gfsu_,gfsv_,cu_,cv_)
        , dof_exchanger(std::make_shared<BorderDOFExchanger>(*this))
        , local_assembler(lop_, cu_, cv_,dof_exchanger)
        , backend(mb_)
      {}
 
      GridOperator(const GFSU & gfsu_, const GFSV & gfsv_, LOP & lop_, const MB& mb_ = MB())
        : global_assembler(gfsu_,gfsv_)
        , dof_exchanger(std::make_shared<BorderDOFExchanger>(*this))
        , local_assembler(lop_,dof_exchanger)
        , backend(mb_)
      {}
 
      const GFSU& trialGridFunctionSpace() const
      {
        return global_assembler.trialGridFunctionSpace();
      }
 
      const GFSV& testGridFunctionSpace() const
      {
        return global_assembler.testGridFunctionSpace();
      }
 
      typename GFSU::Traits::SizeType globalSizeU () const
      {
        return trialGridFunctionSpace().globalSize();
      }
 
      typename GFSV::Traits::SizeType globalSizeV () const
      {
        return testGridFunctionSpace().globalSize();
      }
 
      Assembler & assembler() { return global_assembler; }
 
      const Assembler & assembler() const { return global_assembler; }
 
      LocalAssembler & localAssembler() const { return local_assembler; }
 
 
      template <typename GridOperatorTuple>
      struct SetupGridOperator
      {
        SetupGridOperator()
          : index(0), size(std::tuple_size<GridOperatorTuple>::value) {}
 
        template <typename T>
        void visit(T& elem) {
          elem.localAssembler().preProcessing(index == 0);
          elem.localAssembler().postProcessing(index == size-1);
          ++index;
        }
 
        int index;
        const int size;
      };
 
      template<typename GridOperatorTuple>
      static void setupGridOperators(GridOperatorTuple tuple)
      {
        SetupGridOperator<GridOperatorTuple> setup_visitor;
        Hybrid::forEach(tuple,
                        [&](auto &el) { setup_visitor.visit(el); });
      }
 
      template<typename F, typename X>
      void interpolate (const X& xold, F& f, X& x) const
      {
        // If xold is alias of x it is rewritten by interpolate. Initial guess is then stored in tmp
        std::optional<X> tmp;
        const X& source = (&xold == &x) ? tmp.emplace(xold) : xold;
 
        // Interpolate f into grid function space and set corresponding coefficients
        Dune::PDELab::interpolate(f,global_assembler.trialGridFunctionSpace(),x);
 
        // Copy non-constrained dofs from old time step
        Dune::PDELab::copy_nonconstrained_dofs(local_assembler.trialConstraints(),source,x);
      }
 
      void fill_pattern(Pattern & p) const
      {
        typedef typename LocalAssembler::LocalPatternAssemblerEngine PatternEngine;
        PatternEngine & pattern_engine = local_assembler.localPatternAssemblerEngine(p);
        global_assembler.assemble(pattern_engine);
      }
 
      void residual(const Domain & x, Range & r) const
      {
        typedef typename LocalAssembler::LocalResidualAssemblerEngine ResidualEngine;
        ResidualEngine & residual_engine = local_assembler.localResidualAssemblerEngine(r,x);
        global_assembler.assemble(residual_engine);
      }
 
      void jacobian(const Domain & x, Jacobian & a) const
      {
        typedef typename LocalAssembler::LocalJacobianAssemblerEngine JacobianEngine;
        JacobianEngine & jacobian_engine = local_assembler.localJacobianAssemblerEngine(a,x);
        global_assembler.assemble(jacobian_engine);
      }
 
      void jacobian_apply(const Domain & update, Range & result) const
      {
        if (not local_assembler.localOperator().isLinear)
          DUNE_THROW(Dune::Exception, "Your trying to use a non linear jacobian apply for a linear problem.");
        global_assembler.assemble(local_assembler.localJacobianApplyAssemblerEngine(update, result));
      }
 
      void jacobian_apply(const Domain & solution, const Domain & update, Range & result) const
      {
        if (local_assembler.localOperator().isLinear)
          DUNE_THROW(Dune::Exception, "Your trying to use a non linear jacobian apply for a linear problem.");
        global_assembler.assemble(local_assembler.localJacobianApplyAssemblerEngine(solution, update, result));
      }
 
      [[deprecated("nonlinear_jacobian_apply(x,z,r) is deprecated. Please use jacobian_apply(solution, update, result) instead!")]]
      void nonlinear_jacobian_apply(const Domain & solution, const Domain & update, Range & result) const
      {
        if (local_assembler.localOperator().isLinear)
          DUNE_THROW(Dune::Exception, "Your trying to use a non linear jacobian apply for a linear problem.");
        global_assembler.assemble(local_assembler.localJacobianApplyAssemblerEngine(solution, update, result));
      }
 
      void make_consistent(Jacobian& a) const
      {
        dof_exchanger->accumulateBorderEntries(*this,a);
      }
 
      void update()
      {
        // the DOF exchanger has matrix information, so we need to update it
        dof_exchanger->update(*this);
      }
 
      const typename Traits::MatrixBackend& matrixBackend() const
      {
        return backend;
      }
 
    private:
      Assembler global_assembler;
      std::shared_ptr<BorderDOFExchanger> dof_exchanger;
 
      mutable LocalAssembler local_assembler;
      MB backend;
 
    };
 
  }
}
#endif // DUNE_PDELAB_GRIDOPERATOR_GRIDOPERATOR_HH