residualengine.hh

#ifndef DUNE_PDELAB_GRIDOPERATOR_ONESTEP_RESIDUALENGINE_HH
#define DUNE_PDELAB_GRIDOPERATOR_ONESTEP_RESIDUALENGINE_HH
 
#include <dune/pdelab/gridoperator/onestep/enginebase.hh>
 
#include <cmath>
 
namespace Dune{
  namespace PDELab{
 
    template<typename OSLA>
    class OneStepLocalResidualAssemblerEngine
      : public OneStepLocalAssemblerEngineBase<OSLA,
                                               typename OSLA::LocalAssemblerDT0::LocalResidualAssemblerEngine,
                                               typename OSLA::LocalAssemblerDT1::LocalResidualAssemblerEngine
                                               >
    {
 
      typedef OneStepLocalAssemblerEngineBase<OSLA,
                                              typename OSLA::LocalAssemblerDT0::LocalResidualAssemblerEngine,
                                              typename OSLA::LocalAssemblerDT1::LocalResidualAssemblerEngine
                                              > BaseT;
 
      using BaseT::la;
      using BaseT::lae0;
      using BaseT::lae1;
      using BaseT::implicit;
      using BaseT::setLocalAssemblerEngineDT0;
      using BaseT::setLocalAssemblerEngineDT1;
 
    public:
      typedef OSLA OneStepLocalAssembler;
 
      typedef typename OSLA::LocalAssemblerDT0 LocalAssemblerDT0;
      typedef typename OSLA::LocalAssemblerDT1 LocalAssemblerDT1;
 
      typedef typename LocalAssemblerDT0::LocalResidualAssemblerEngine ResidualEngineDT0;
      typedef typename LocalAssemblerDT1::LocalResidualAssemblerEngine ResidualEngineDT1;
 
      typedef typename OSLA::Traits::Residual Residual;
 
      typedef typename OSLA::Traits::Solution Solution;
 
      typedef typename OSLA::Real Real;
 
      typedef OSLA LocalAssembler;
 
      OneStepLocalResidualAssemblerEngine(const LocalAssembler & local_assembler_)
        : BaseT(local_assembler_)
        , invalid_residual(nullptr)
        , invalid_solution(nullptr)
        , residual_0(invalid_residual)
        , residual_1(invalid_residual)
        , const_residual_0(invalid_residual)
        , const_residual_1(invalid_residual)
        , solution(invalid_solution)
      {}
 
      void setSolution(const Solution & solution_)
      {
        solution = &solution_;
      }
 
      void setConstResidual(const Residual &const_residual_)
      {
        const_residual_0 = &const_residual_;
        const_residual_1 = &const_residual_;
      }
 
      void setResidual(Residual & residual_)
      {
        residual_0 = &residual_;
        residual_1 = &residual_;
 
        // Initialize the engines of the two wrapped local assemblers
        assert(solution != invalid_solution);
        setLocalAssemblerEngineDT0(la.la0.localResidualAssemblerEngine(*residual_0,*solution));
        setLocalAssemblerEngineDT1(la.la1.localResidualAssemblerEngine(*residual_1,*solution));
      }
 
      void setConstResiduals(const Residual &const_residual_0_, const Residual &const_residual_1_)
      {
        const_residual_0 = &const_residual_0_;
        const_residual_1 = &const_residual_1_;
      }
 
      void setResiduals(Residual & residual_0_, Residual & residual_1_)
      {
        residual_0 = &residual_0_;
        residual_1 = &residual_1_;
 
        // Initialize the engines of the two wrapped local assemblers
        assert(solution != invalid_solution);
        setLocalAssemblerEngineDT0(la.la0.localResidualAssemblerEngine(*residual_0,*solution));
        setLocalAssemblerEngineDT1(la.la1.localResidualAssemblerEngine(*residual_1,*solution));
      }
 
      void setWeights()
      {
        la.la0.setWeight(b_rr * la.dt_factor0);
        la.la1.setWeight(la.dt_factor1);
      }
 
      void preAssembly()
      {
        lae0->preAssembly();
        lae1->preAssembly();
 
        // Extract the coefficients of the time step scheme
        b_rr = la.osp_method->b(la.stage,la.stage);
        d_r = la.osp_method->d(la.stage);
        using std::abs;
        implicit = abs(b_rr) > 1e-6;
 
        // prepare local operators for stage
        la.la0.setTime(la.time + d_r * la.dt);
        la.la1.setTime(la.time + d_r * la.dt);
 
        setWeights();
      }
 
      template<typename GFSU, typename GFSV>
      void postAssembly(const GFSU& gfsu, const GFSV& gfsv)
      {
 
        // Update residual vectors with constant part
        assert(const_residual_0 != invalid_residual);
        assert(const_residual_1 != invalid_residual);
        *residual_0 += *const_residual_0;
        if(residual_0 != residual_1){
          assert(const_residual_0 != const_residual_1);
          *residual_1 += *const_residual_1;
        }
 
        lae0->postAssembly(gfsu,gfsv);
        lae1->postAssembly(gfsu,gfsv);
      }
 
 
    private:
 
      Residual * const invalid_residual;
 
      Solution * const invalid_solution;
 
      Residual * residual_0;
      Residual * residual_1;
 
      const Residual * const_residual_0;
      const Residual * const_residual_1;
 
      const Solution * solution;
 
      Real b_rr, d_r;
 
    }; // End of class OneStepLocalResidualAssemblerEngine
 
  }
}
 
#endif // DUNE_PDELAB_GRIDOPERATOR_ONESTEP_RESIDUALENGINE_HH