tmixed_upscaledc.hh

#pragma once

#include<dune/pdelab/common/quadraturerules.hh>
#include<dune/pdelab/finiteelement/localbasiscache.hh>
#include<dune/pdelab/localoperator/defaultimp.hh>
#include<dune/pdelab/localoperator/pattern.hh>
#include<dune/pdelab/localoperator/flags.hh>
#include<dune/pdelab/localoperator/variablefactories.hh>
#include<dune/pdelab/gridfunctionspace/gridfunctionspace.hh>





template<typename Param, typename PAdaptor, typename CGFS, typename CContainer, typename RF>
    class TMIXED_UpscaledC :
      public Dune::PDELab::FullVolumePattern,
      public  Dune::PDELab::FullSkeletonPattern,
      public Dune::PDELab::LocalOperatorDefaultFlags,
      public Dune::PDELab::NumericalJacobianVolume<TMIXED_UpscaledC<Param, PAdaptor, CGFS, CContainer, RF> >,
      public Dune::PDELab::NumericalJacobianSkeleton<TMIXED_UpscaledC<Param, PAdaptor, CGFS, CContainer, RF> >
    {
    private:
      Param& param; // parameter functions
      PAdaptor& pAdapt;

      //typedef typename CGFS::Traits::Type FEM;
      typedef typename CGFS::Traits::FiniteElementType::Traits::LocalBasisType LocalBasis;
      Dune::PDELab::LocalBasisCache<LocalBasis> cache;

      typedef Dune::PDELab::LocalFunctionSpace<CGFS> CLFS;
      mutable CLFS cLfs;
      typedef Dune::PDELab::LFSIndexCache<CLFS> CLFSCache;
      mutable CLFSCache cLfsCache;
      typedef typename CContainer::template ConstLocalView<CLFSCache> CView;
      mutable CView cView;
      mutable std::vector<typename LocalBasis::Traits::RangeFieldType> cOldLocal;

      //FOR DETECTING ENTITY changes
      typedef typename CGFS::Traits::GridViewType::template Codim<0>::Entity   ENT;
      mutable ENT OldEnt;

    public:
      // pattern assembly flags
      enum { doPatternVolume = true };
      enum { doPatternSkeleton = true};

      // residual assembly flags
      enum { doAlphaVolume = true };
      //enum { doLambdaVolume = true };
      enum { doAlphaSkeleton = true};


      TMIXED_UpscaledC (Param& param_, PAdaptor& pAdapt_, CGFS& cGfs_, CContainer& cVOld_) :
        param(param_), pAdapt(pAdapt_),
        cLfs(cGfs_), cLfsCache(cLfs), cView(cVOld_), cOldLocal(cGfs_.maxLocalSize())
      {
      }


      template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
      void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
      {
        std::cout << "hi alpha" << std::endl;
        auto c = x(lfsu,0);
        if(eg.entity() != OldEnt)
        {
          OldEnt = eg.entity();
          cLfs.bind(eg.entity());
          cLfsCache.update();
          cView.bind(cLfsCache);
          cView.read(cOldLocal);
          cView.unbind();
        }
        auto cOld = cOldLocal[lfsu.localIndex(0)];

        r.accumulate(lfsv,0, (c-cOld)/param.time.dt*eg.geometry().volume());
        std::cout << "bye alpha" << std::endl;
      } //alpha volume


      //! residual contribution from skeleton terms
      template<typename IG, typename LFSU, typename X,
               typename LFSV, typename R>
      void alpha_skeleton (const IG& ig,
             const LFSU& lfsu_i, const X& x_i, const LFSV& lfsv_i,
             const LFSU& lfsu_o, const X& x_o, const LFSV& lfsv_o,
             R& r_i, R& r_o) const
      {
        std::cout << "hi skeleton" << std::endl;
        // inside and outside cells
        auto cell_inside = ig.inside();
        auto cell_outside = ig.outside();

        // inside and outside cell geometries
        auto insidegeo = cell_inside.geometry();
        auto outsidegeo = cell_outside.geometry();

        // cell centers in global coordinates
        auto inside_global = insidegeo.center();
        auto outside_global = outsidegeo.center();

        // distance between the two cell centers
        inside_global -= outside_global;
        auto distance = inside_global.two_norm();

        // face volume for integration
        //auto facegeo = ig.geometry();
        //auto face_volume = facegeo.volume();

        // contribution to residual on inside and outside elements
        auto dudn = (x_o(lfsu_i,0)-x_i(lfsu_o,0))/distance;
        auto DiffusionCoeff = 0.5*(pAdapt.evalKf(inside_global) + pAdapt.evalKf(outside_global));
        r_i.accumulate(lfsv_i,0,-dudn*DiffusionCoeff);
        r_o.accumulate(lfsv_o,0, dudn*DiffusionCoeff);
        std::cout << "bye skeleton" << std::endl;
      }


    };