jacobianapplyengine.hh

#ifndef DUNE_PDELAB_GRIDOPERATOR_DEFAULT_JACOBIANAPPLYENGINE_HH
#define DUNE_PDELAB_GRIDOPERATOR_DEFAULT_JACOBIANAPPLYENGINE_HH
 
#include <dune/pdelab/gridfunctionspace/localvector.hh>
#include <dune/pdelab/gridoperator/common/assemblerutilities.hh>
#include <dune/pdelab/gridoperator/common/localassemblerenginebase.hh>
#include <dune/pdelab/constraints/common/constraints.hh>
#include <dune/pdelab/localoperator/callswitch.hh>
 
namespace Dune{
  namespace PDELab{
 
    template<typename LA>
    class DefaultLocalJacobianApplyAssemblerEngine
      : public LocalAssemblerEngineBase
    {
    public:
 
      template<typename TrialConstraintsContainer, typename TestConstraintsContainer>
      bool needsConstraintsCaching(const TrialConstraintsContainer& cu, const TestConstraintsContainer& cv) const
      {
        return false;
      }
 
      typedef LA LocalAssembler;
 
      typedef typename LA::LocalOperator LOP;
 
      static constexpr bool isLinear = LOP::isLinear;
 
      typedef typename LA::Traits::Range Range;
      typedef typename Range::ElementType RangeElement;
 
      typedef typename LA::Traits::Domain Domain;
      typedef typename Domain::ElementType DomainElement;
 
      typedef typename LA::LFSU LFSU;
      typedef typename LA::LFSUCache LFSUCache;
      typedef typename LFSU::Traits::GridFunctionSpace GFSU;
      typedef typename LA::LFSV LFSV;
      typedef typename LA::LFSVCache LFSVCache;
      typedef typename LFSV::Traits::GridFunctionSpace GFSV;
 
      typedef typename Domain::template ConstLocalView<LFSUCache> DomainView;
      typedef typename Range::template LocalView<LFSVCache> RangeView;
 
      DefaultLocalJacobianApplyAssemblerEngine(const LocalAssembler & local_assembler_)
        : local_assembler(local_assembler_),
          lop(local_assembler_.localOperator()),
          result_view_inside(local_result_inside,1.0),
          result_view_outside(local_result_outside,1.0)
      {}
 
      bool requireSkeleton() const
      { return local_assembler.doAlphaSkeleton(); }
      bool requireSkeletonTwoSided() const
      { return local_assembler.doSkeletonTwoSided(); }
      bool requireUVVolume() const
      { return local_assembler.doAlphaVolume(); }
      bool requireUVSkeleton() const
      { return local_assembler.doAlphaSkeleton(); }
      bool requireUVBoundary() const
      { return local_assembler.doAlphaBoundary(); }
      bool requireUVVolumePostSkeleton() const
      { return local_assembler.doAlphaVolumePostSkeleton(); }
 
      const LocalAssembler & localAssembler() const
      {
        return local_assembler;
      }
 
      const typename LocalAssembler::Traits::TrialGridFunctionSpaceConstraints& trialConstraints() const
      {
        return localAssembler().trialConstraints();
      }
 
      const typename LocalAssembler::Traits::TestGridFunctionSpaceConstraints& testConstraints() const
      {
        return localAssembler().testConstraints();
      }
 
      void setSolution(const Domain & solution_)
      {
        global_solution_view_inside.attach(solution_);
        global_solution_view_outside.attach(solution_);
      }
 
      void setUpdate(const Domain & update_)
      {
        global_update_view_inside.attach(update_);
        global_update_view_outside.attach(update_);
      }
 
      void setResult(Range & result_)
      {
        global_result_view_inside.attach(result_);
        global_result_view_outside.attach(result_);
      }
 
      template<typename EG, typename LFSUC, typename LFSVC>
      void onBindLFSUV(const EG & eg, const LFSUC & lfsu_cache, const LFSVC & lfsv_cache)
      {
        global_solution_view_inside.bind(lfsu_cache);
        // In the linear case the solution is not needed and we drop this part
        if (not isLinear)
          local_solution_inside.resize(lfsu_cache.size());
        global_update_view_inside.bind(lfsu_cache);
        local_update_inside.resize(lfsu_cache.size());
      }
 
      template<typename EG, typename LFSVC>
      void onBindLFSV(const EG & eg, const LFSVC & lfsv_cache)
      {
        global_result_view_inside.bind(lfsv_cache);
        local_result_inside.assign(lfsv_cache.size(),0.0);
      }
 
      template<typename IG, typename LFSUC, typename LFSVC>
      void onBindLFSUVInside(const IG & ig, const LFSUC & lfsu_cache, const LFSVC & lfsv_cache)
      {
        global_solution_view_inside.bind(lfsu_cache);
        // In the linear case the solution is not needed and we drop this part
        if (not isLinear)
          local_solution_inside.resize(lfsu_cache.size());
        global_update_view_inside.bind(lfsu_cache);
        local_update_inside.resize(lfsu_cache.size());
      }
 
      template<typename IG, typename LFSUC, typename LFSVC>
      void onBindLFSUVOutside(const IG & ig,
                              const LFSUC & lfsu_s_cache, const LFSVC & lfsv_s_cache,
                              const LFSUC & lfsu_n_cache, const LFSVC & lfsv_n_cache)
      {
        global_solution_view_outside.bind(lfsu_n_cache);
        // In the linear case the solution is not needed and we drop this part
        if (not isLinear)
          local_solution_outside.resize(lfsu_n_cache.size());
        global_update_view_outside.bind(lfsu_n_cache);
        local_update_outside.resize(lfsu_n_cache.size());
      }
 
      template<typename IG, typename LFSVC>
      void onBindLFSVInside(const IG & ig, const LFSVC & lfsv_cache)
      {
        global_result_view_inside.bind(lfsv_cache);
        local_result_inside.assign(lfsv_cache.size(),0.0);
      }
 
      template<typename IG, typename LFSVC>
      void onBindLFSVOutside(const IG & ig,
                             const LFSVC & lfsv_s_cache,
                             const LFSVC & lfsv_n_cache)
      {
        global_result_view_outside.bind(lfsv_n_cache);
        local_result_outside.assign(lfsv_n_cache.size(),0.0);
      }
 
 
      template<typename EG, typename LFSVC>
      void onUnbindLFSV(const EG & eg, const LFSVC & lfsv_cache)
      {
        global_result_view_inside.add(local_result_inside);
        global_result_view_inside.commit();
      }
 
      template<typename IG, typename LFSVC>
      void onUnbindLFSVInside(const IG & ig, const LFSVC & lfsv_cache)
      {
        global_result_view_inside.add(local_result_inside);
        global_result_view_inside.commit();
      }
 
      template<typename IG, typename LFSVC>
      void onUnbindLFSVOutside(const IG & ig,
                               const LFSVC & lfsv_s_cache,
                               const LFSVC & lfsv_n_cache)
      {
        global_result_view_outside.add(local_result_outside);
        global_result_view_outside.commit();
      }
 
      template<typename LFSUC>
      void loadCoefficientsLFSUInside(const LFSUC & lfsu_s_cache)
      {
        // In the linear case the solution is not needed and we drop this part
        if (not isLinear)
          global_solution_view_inside.read(local_solution_inside);
        global_update_view_inside.read(local_update_inside);
      }
      template<typename LFSUC>
      void loadCoefficientsLFSUOutside(const LFSUC & lfsu_n_cache)
      {
        // In the linear case the solution is not needed and we drop this part
        if (not isLinear)
          global_solution_view_outside.read(local_solution_outside);
        global_update_view_outside.read(local_update_outside);
      }
      template<typename LFSUC>
      void loadCoefficientsLFSUCoupling(const LFSUC & lfsu_c_cache)
      {
        DUNE_THROW(Dune::NotImplemented,"No coupling lfsu available for ");
      }
 
 
      void postAssembly(const GFSU& gfsu, const GFSV& gfsv)
      {
        if(local_assembler.doPostProcessing())
          Dune::PDELab::constrain_residual(local_assembler.testConstraints(),
                                           global_result_view_inside.container());
      }
 
 
 
      template <typename EG> bool skipEntity(const EG &eg)
      {
        return localAssembler().skipEntity(eg) ||
               (LocalAssembler::isNonOverlapping &&
                eg.entity().partitionType() != Dune::InteriorEntity);
      }
 
      template<typename IG>
      bool skipIntersection(const IG & ig)
      {
        return localAssembler().skipIntersection(ig);
      }
 
      template<typename EG, typename LFSUC, typename LFSVC>
      void assembleUVVolume(const EG & eg, const LFSUC & lfsu_cache, const LFSVC & lfsv_cache)
      {
        result_view_inside.setWeight(local_assembler.weight());
        Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaVolume>::
          jacobian_apply_volume(lop,eg,lfsu_cache.localFunctionSpace(),local_solution_inside,local_update_inside,lfsv_cache.localFunctionSpace(),result_view_inside);
      }
 
      template<typename IG, typename LFSUC, typename LFSVC>
      void assembleUVSkeleton(const IG & ig, const LFSUC & lfsu_s_cache, const LFSVC & lfsv_s_cache,
                              const LFSUC & lfsu_n_cache, const LFSVC & lfsv_n_cache)
      {
        result_view_inside.setWeight(local_assembler.weight());
        result_view_outside.setWeight(local_assembler.weight());
        Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaSkeleton>::
          jacobian_apply_skeleton(lop,ig,
                                  lfsu_s_cache.localFunctionSpace(),local_solution_inside,local_update_inside,lfsv_s_cache.localFunctionSpace(),
                                  lfsu_n_cache.localFunctionSpace(),local_solution_outside,local_update_outside,lfsv_n_cache.localFunctionSpace(),
                                  result_view_inside,result_view_outside);
      }
 
      template<typename IG, typename LFSUC, typename LFSVC>
      void assembleUVBoundary(const IG & ig, const LFSUC & lfsu_s_cache, const LFSVC & lfsv_s_cache)
      {
        result_view_inside.setWeight(local_assembler.weight());
        Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaBoundary>::
          jacobian_apply_boundary(lop,ig,lfsu_s_cache.localFunctionSpace(),local_solution_inside,local_update_inside,lfsv_s_cache.localFunctionSpace(),result_view_inside);
      }
 
      template<typename IG, typename LFSUC, typename LFSVC>
      static void assembleUVEnrichedCoupling(const IG & ig,
                                             const LFSUC & lfsu_s_cache, const LFSVC & lfsv_s_cache,
                                             const LFSUC & lfsu_n_cache, const LFSVC & lfsv_n_cache,
                                             const LFSUC & lfsu_coupling_cache, const LFSVC & lfsv_coupling_cache)
      {
        DUNE_THROW(Dune::NotImplemented,"Assembling of coupling spaces is not implemented for ");
      }
 
      template<typename EG, typename LFSUC, typename LFSVC>
      void assembleUVVolumePostSkeleton(const EG & eg, const LFSUC & lfsu_cache, const LFSVC & lfsv_cache)
      {
        result_view_inside.setWeight(local_assembler.weight());
        Dune::PDELab::LocalAssemblerCallSwitch<LOP,LOP::doAlphaVolumePostSkeleton>::
          jacobian_apply_volume_post_skeleton(lop,eg,lfsu_cache.localFunctionSpace(),local_solution_inside,local_update_inside,lfsv_cache.localFunctionSpace(),result_view_inside);
      }
 
 
    private:
      const LocalAssembler & local_assembler;
 
      const LOP & lop;
 
      DomainView global_solution_view_inside;
      DomainView global_solution_view_outside;
 
      DomainView global_update_view_inside;
      DomainView global_update_view_outside;
 
      RangeView global_result_view_inside;
      RangeView global_result_view_outside;
 
      typedef Dune::PDELab::TrialSpaceTag LocalTrialSpaceTag;
      typedef Dune::PDELab::TestSpaceTag LocalTestSpaceTag;
 
      typedef Dune::PDELab::LocalVector<DomainElement, LocalTrialSpaceTag> DomainVector;
      typedef Dune::PDELab::LocalVector<RangeElement, LocalTestSpaceTag> RangeVector;
 
      DomainVector local_solution_inside;
      DomainVector local_solution_outside;
      DomainVector local_update_inside;
      DomainVector local_update_outside;
      RangeVector local_result_inside;
      RangeVector local_result_outside;
      typename RangeVector::WeightedAccumulationView result_view_inside;
      typename RangeVector::WeightedAccumulationView result_view_outside;
 
    }; // End of class DefaultLocalJacobianAssemblerEngine
 
  }
}
#endif // DUNE_PDELAB_GRIDOPERATOR_DEFAULT_JACOBIANAPPLYENGINE_HH