jacobianengine.hh

#ifndef DUNE_PDELAB_GRIDOPERATOR_ONESTEP_JACOBIANENGINE_HH
#define DUNE_PDELAB_GRIDOPERATOR_ONESTEP_JACOBIANENGINE_HH
 
#include <dune/pdelab/gridoperator/onestep/enginebase.hh>
#include <cmath>
 
namespace Dune{
  namespace PDELab{
 
    template<typename OSLA>
    class OneStepLocalJacobianAssemblerEngine
      : public OneStepLocalAssemblerEngineBase<OSLA,
                                               typename OSLA::LocalAssemblerDT0::LocalJacobianAssemblerEngine,
                                               typename OSLA::LocalAssemblerDT1::LocalJacobianAssemblerEngine
                                               >
    {
 
      typedef OneStepLocalAssemblerEngineBase<OSLA,
                                              typename OSLA::LocalAssemblerDT0::LocalJacobianAssemblerEngine,
                                              typename OSLA::LocalAssemblerDT1::LocalJacobianAssemblerEngine
                                              > BaseT;
 
      using BaseT::la;
      using BaseT::lae0;
      using BaseT::lae1;
      using BaseT::implicit;
      using BaseT::setLocalAssemblerEngineDT0;
      using BaseT::setLocalAssemblerEngineDT1;
    public:
      typedef OSLA LocalAssembler;
 
      typedef typename OSLA::LocalAssemblerDT0 LocalAssemblerDT0;
      typedef typename OSLA::LocalAssemblerDT1 LocalAssemblerDT1;
 
      typedef typename LocalAssemblerDT0::LocalJacobianAssemblerEngine JacobianEngineDT0;
      typedef typename LocalAssemblerDT1::LocalJacobianAssemblerEngine JacobianEngineDT1;
 
      typedef typename OSLA::Traits::Jacobian Jacobian;
 
      typedef typename OSLA::Traits::Solution Solution;
 
      typedef typename OSLA::Real Real;
 
      OneStepLocalJacobianAssemblerEngine(const LocalAssembler & local_assembler_)
        : BaseT(local_assembler_),
          invalid_jacobian(nullptr),
          invalid_solution(nullptr),
          jacobian(invalid_jacobian), solution(invalid_solution)
      {}
 
 
      void setSolution(const Solution & solution_)
      {
        solution = &solution_;
      }
 
      void setJacobian(Jacobian & jacobian_)
      {
        jacobian = &jacobian_;
 
        assert(solution != invalid_solution);
 
        // Initialize the engines of the two wrapped local assemblers
        setLocalAssemblerEngineDT0(la.la0.localJacobianAssemblerEngine(*jacobian,*solution));
        setLocalAssemblerEngineDT1(la.la1.localJacobianAssemblerEngine(*jacobian,*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);
 
        // Here we only want to know whether this stage is implicit
        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)
      {
        lae0->postAssembly(gfsu,gfsv);
        lae1->postAssembly(gfsu,gfsv);
      }
 
    private:
 
      Jacobian * const invalid_jacobian;
 
      Solution * const invalid_solution;
 
      Jacobian * jacobian;
 
      const Solution * solution;
 
      Real b_rr, d_r;
 
    }; // End of class OneStepLocalJacobianAssemblerEngine
 
  }
}
 
#endif // DUNE_PDELAB_GRIDOPERATOR_ONESTEP_JACOBIANENGINE_HH