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 xg = eg.geometry().center();
        auto cOld = cOldLocal[lfsu.localIndex(0)];
        auto phic = pAdapt.evalPhiC(xg);
        auto phicOld = pAdapt.evalPhiCOld(xg);
        auto phiSolid = pAdapt.evalPhiSolid(xg);
        auto phiSolidOld = pAdapt.evalPhiSolidOld(xg);
        r.accumulate(lfsv,0, eg.geometry().volume()*
          ((phic*c - phicOld*cOld + param.phys.uAst*(phiSolid-phiSolidOld))/param.time.dt));
        //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();
        auto facegeo = ig.geometry();
        auto face_global = facegeo.center();

        //std::cout << "face: " << face_global[0] << "inside:" <<
        //    inside_global[0] << "outside: " << outside_global[0] << std::endl;


        // distance between the two cell centers
        auto distanceVector = inside_global;
        distanceVector -= outside_global;
        auto distance = distanceVector.two_norm();
        bool insideIsUpwind = (inside_global[0] < outside_global[0]);


        auto insideKc = pAdapt.evalKc(inside_global);
        auto outsideKc = pAdapt.evalKc(outside_global);
        auto insideDxP = pAdapt.evalDxP(inside_global);
        auto outsideDxP = pAdapt.evalDxP(outside_global);
        auto insideC = x_i(lfsu_i,0);
        auto outsideC = x_o(lfsu_o,0);

        if(distance > 0.5)
        {
          insideIsUpwind = !insideIsUpwind;
          distance = 1-distance;
          //std::cout << "Edge detected" << distance << " "
          //  <<insideC<< " "<<outsideC<< " " << insideIsUpwind << std::endl;
        }
        //std::cout << inside_global[0] << " " << outside_global[0] << " " << insideIsUpwind << std::endl;

        //auto flux = outsideKc*outsideDxP*outsideC - insideKc*insideDxP*insideC;
        auto flux = outsideKc * outsideDxP * outsideC;
        if(insideIsUpwind)
        {
          flux = - insideKc*insideDxP*insideC;
        }

        // face volume for integratio
        //auto face_volume = facegeo.volume();

        // contribution to residual on inside and outside elements
        auto dudn = (outsideC-insideC)/distance;
        RF peclet = 5;  //TODO THIS HAS TO BE FIXED IN OPTIONS!!!!!!!!!!
        auto DiffusionCoeff = 1/peclet * 0.5*(pAdapt.evalPhiC(inside_global) + pAdapt.evalPhiC(outside_global));
        r_i.accumulate(lfsv_i,0,-dudn*DiffusionCoeff + flux);
        r_o.accumulate(lfsv_o,0, dudn*DiffusionCoeff - flux);
        //std::cout << "bye skeleton" << std::endl;
      }


    };