localassembler.hh

// -*- tab-width: 2; indent-tabs-mode: nil -*-
// vi: set et ts=2 sw=2 sts=2:
#ifndef DUNE_PDELAB_GRIDOPERATOR_FASTDG_LOCALASSEMBLER_HH
#define DUNE_PDELAB_GRIDOPERATOR_FASTDG_LOCALASSEMBLER_HH
 
#include <memory>
 
#include <dune/common/shared_ptr.hh>
 
#include <dune/typetree/typetree.hh>
 
#include <dune/pdelab/gridoperator/fastdg/residualengine.hh>
#include <dune/pdelab/gridoperator/fastdg/patternengine.hh>
#include <dune/pdelab/gridoperator/fastdg/jacobianengine.hh>
#include <dune/pdelab/gridoperator/fastdg/jacobianapplyengine.hh>
#include <dune/pdelab/gridoperator/common/assemblerutilities.hh>
#include <dune/pdelab/gridfunctionspace/lfsindexcache.hh>
 
namespace Dune{
  namespace PDELab{
 
    template<typename GO, typename LOP, bool nonoverlapping_mode = false>
    class FastDGLocalAssembler :
      public Dune::PDELab::LocalAssemblerBase<typename GO::Traits::MatrixBackend,
                                              typename GO::Traits::TrialGridFunctionSpaceConstraints,
                                              typename GO::Traits::TestGridFunctionSpaceConstraints>
    {
    public:
 
      typedef Dune::PDELab::LocalAssemblerTraits<GO> Traits;
 
      typedef typename Traits::Residual::ElementType RangeField;
      typedef RangeField Real;
 
      typedef typename Traits::TrialGridFunctionSpace GFSU;
      typedef typename Traits::TestGridFunctionSpace GFSV;
 
      typedef typename Traits::TrialGridFunctionSpaceConstraints CU;
      typedef typename Traits::TestGridFunctionSpaceConstraints CV;
 
      typedef Dune::PDELab::LocalAssemblerBase<typename Traits::MatrixBackend,CU,CV> Base;
 
      typedef LOP LocalOperator;
 
      static const bool isNonOverlapping = nonoverlapping_mode;
 
      // Types of local function spaces
      typedef Dune::PDELab::LocalFunctionSpace<GFSU, Dune::PDELab::TrialSpaceTag> LFSU;
      typedef Dune::PDELab::LocalFunctionSpace<GFSV, Dune::PDELab::TestSpaceTag> LFSV;
      typedef LFSIndexCache<LFSU,CU,true> LFSUCache;
      typedef LFSIndexCache<LFSV,CV,true> LFSVCache;
 
 
      typedef FastDGLocalPatternAssemblerEngine<FastDGLocalAssembler> LocalPatternAssemblerEngine;
      typedef FastDGLocalResidualAssemblerEngine<FastDGLocalAssembler> LocalResidualAssemblerEngine;
      typedef FastDGLocalJacobianAssemblerEngine<FastDGLocalAssembler> LocalJacobianAssemblerEngine;
      typedef FastDGLocalJacobianApplyAssemblerEngine<FastDGLocalAssembler> LocalJacobianApplyAssemblerEngine;
 
      // friend declarations such that engines are able to call scatter_jacobian() and add_entry() from base class
      friend class FastDGLocalPatternAssemblerEngine<FastDGLocalAssembler>;
 
      FastDGLocalAssembler (LOP & lop_, std::shared_ptr<typename GO::BorderDOFExchanger> border_dof_exchanger)
        : lop(stackobject_to_shared_ptr(lop_)),
          weight_(1.0),
          doPreProcessing_(true),
          doPostProcessing_(true),
          pattern_engine(*this,border_dof_exchanger), residual_engine(*this), jacobian_engine(*this)
        , jacobian_apply_engine(*this)
        , _reconstruct_border_entries(isNonOverlapping)
      {}
 
      FastDGLocalAssembler (LOP & lop_, const CU& cu_, const CV& cv_,
                            std::shared_ptr<typename GO::BorderDOFExchanger> border_dof_exchanger)
        : Base(cu_, cv_),
          lop(stackobject_to_shared_ptr(lop_)),
          weight_(1.0),
          doPreProcessing_(true),
          doPostProcessing_(true),
          pattern_engine(*this,border_dof_exchanger), residual_engine(*this), jacobian_engine(*this)
        , jacobian_apply_engine(*this)
        , _reconstruct_border_entries(isNonOverlapping)
      {}
 
      LOP &localOperator()
      {
        return *lop;
      }
 
      const LOP &localOperator() const
      {
        return *lop;
      }
 
      template<class EG>
      bool skipEntity(const EG & eg) const
      {
        static_assert(
            models<Impl::HasDoSkipEntity, LOP>(),
            "Your local operator does not provide the 'doSkipEntity' flag. "
            "If you are porting a previous implementation, set the flag to "
            "false");
        return Dune::PDELab::LocalAssemblerCallSwitch<
            LOP, LOP::doSkipEntity>::skip_entity(*lop, eg);
      }
 
      template<class IG>
      bool skipIntersection(const IG & ig) const
      {
        static_assert(models<Impl::HasDoSkipIntersection, LOP>(),
                      "Your local operator does not provide the "
                      "'doSkipIntersection' flag. "
                      "If you are porting a previous implementation, set the "
                      "flag to false");
        return Dune::PDELab::LocalAssemblerCallSwitch<
            LOP, LOP::doSkipIntersection>::skip_intersection(*lop, ig);
      }
 
      void setTime(Real time_)
      {
        lop->setTime(time_);
      }
 
      RangeField weight() const
      {
        return weight_;
      }
 
      void setWeight(RangeField weight){
        weight_ = weight;
      }
 
      void preStage (Real time_, int r_) { lop->preStage(time_,r_); }
      void preStep (Real time_, Real dt_, std::size_t stages_){ lop->preStep(time_,dt_,stages_); }
      void postStep (){ lop->postStep(); }
      void postStage (){ lop->postStage(); }
      Real suggestTimestep (Real dt) const{return lop->suggestTimestep(dt); }
 
      bool reconstructBorderEntries() const
      {
        return _reconstruct_border_entries;
      }
 
 
      LocalPatternAssemblerEngine & localPatternAssemblerEngine
      (typename Traits::MatrixPattern & p)
      {
        pattern_engine.setPattern(p);
        return pattern_engine;
      }
 
      LocalResidualAssemblerEngine & localResidualAssemblerEngine
      (typename Traits::Residual & r, const typename Traits::Solution & x)
      {
        residual_engine.setResidual(r);
        residual_engine.setSolution(x);
        return residual_engine;
      }
 
      LocalJacobianAssemblerEngine & localJacobianAssemblerEngine
      (typename Traits::Jacobian & a, const typename Traits::Solution & x)
      {
        jacobian_engine.setJacobian(a);
        jacobian_engine.setSolution(x);
        return jacobian_engine;
      }
 
      LocalJacobianApplyAssemblerEngine & localJacobianApplyAssemblerEngine
      (const typename Traits::Domain & update, typename Traits::Range & result)
      {
        jacobian_apply_engine.setUpdate(update);
        jacobian_apply_engine.setResult(result);
        return jacobian_apply_engine;
      }
 
      LocalJacobianApplyAssemblerEngine & localJacobianApplyAssemblerEngine
      (const typename Traits::Domain & solution, const typename Traits::Domain & update, typename Traits::Range & result)
      {
        jacobian_apply_engine.setSolution(solution);
        jacobian_apply_engine.setUpdate(update);
        jacobian_apply_engine.setResult(result);
        return jacobian_apply_engine;
      }
 
 
      static constexpr bool doAlphaVolume() { return LOP::doAlphaVolume; }
      static constexpr bool doLambdaVolume() { return LOP::doLambdaVolume; }
      static constexpr bool doAlphaSkeleton() { return LOP::doAlphaSkeleton; }
      static constexpr bool doLambdaSkeleton() { return LOP::doLambdaSkeleton; }
      static constexpr bool doAlphaBoundary()  { return LOP::doAlphaBoundary; }
      static constexpr bool doLambdaBoundary() { return LOP::doLambdaBoundary; }
      static constexpr bool doAlphaVolumePostSkeleton()  { return LOP::doAlphaVolumePostSkeleton; }
      static constexpr bool doLambdaVolumePostSkeleton() { return LOP::doLambdaVolumePostSkeleton; }
      static constexpr bool doSkeletonTwoSided()  { return LOP::doSkeletonTwoSided; }
      static constexpr bool doPatternVolume()  { return LOP::doPatternVolume; }
      static constexpr bool doPatternSkeleton()  { return LOP::doPatternSkeleton; }
      static constexpr bool doPatternBoundary()  { return LOP::doPatternBoundary; }
      static constexpr bool doPatternVolumePostSkeleton()  { return LOP::doPatternVolumePostSkeleton; }
      static constexpr bool isLinear() { return LOP::isLinear;}
 
 
      bool doPreProcessing() const { return doPreProcessing_; }
 
      void preProcessing(bool v)
      {
        doPreProcessing_ = v;
      }
 
 
      bool doPostProcessing() const { return doPostProcessing_; }
 
      void postProcessing(bool v)
      {
        doPostProcessing_ = v;
      }
 
      template<typename GFSV, typename GC, typename C>
      void set_trivial_rows(const GFSV& gfsv, GC& globalcontainer, const C& c) const
      {
        typedef typename C::const_iterator global_row_iterator;
        for (global_row_iterator cit = c.begin(); cit != c.end(); ++cit)
          globalcontainer.clear_row_block(cit->first,1);
      }
 
      template<typename GFSV, typename GC>
      void set_trivial_rows(const GFSV& gfsv, GC& globalcontainer, const EmptyTransformation& c) const
      {
      }
 
      template<typename GFSV, typename GC>
      void handle_dirichlet_constraints(const GFSV& gfsv, GC& globalcontainer) const
      {
        globalcontainer.flush();
        set_trivial_rows(gfsv,globalcontainer,*this->pconstraintsv);
        globalcontainer.finalize();
      }
 
    private:
 
      std::shared_ptr<LOP> lop;
 
      RangeField weight_;
 
      bool doPreProcessing_;
 
      bool doPostProcessing_;
 
      LocalPatternAssemblerEngine  pattern_engine;
      LocalResidualAssemblerEngine residual_engine;
      LocalJacobianAssemblerEngine jacobian_engine;
      LocalJacobianApplyAssemblerEngine jacobian_apply_engine;
 
      bool _reconstruct_border_entries;
    };
 
  } // end namespace PDELab
} // end namespace Dune
#endif