fem_1p_navierstokes.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>

#include"fem_base_porenetwork.hh"




template<typename Param, typename CHGFS, typename CHContainer, typename NSGFS, typename NSContainer>
    class FEM_1P_NavierStokes :
      public Dune::PDELab::FullVolumePattern,
      public Dune::PDELab::LocalOperatorDefaultFlags,
      public Dune::PDELab::NumericalJacobianVolume<FEM_1P_NavierStokes<Param, CHGFS, CHContainer, NSGFS, NSContainer> >,
      public Dune::PDELab::NumericalJacobianApplyVolume<FEM_1P_NavierStokes<Param, CHGFS, CHContainer, NSGFS, NSContainer> >
    {
    private:
      typedef typename CHContainer::field_type RF;
      typedef typename NSGFS::Traits::GridViewType::template Codim<0>::Entity   ENTITY;

      CHLocalBasisCache<CHGFS> chCache;
      Micro_VXLocalBasisCache<NSGFS> vCache;

      typedef GFS_LocalViewer<CHGFS, CHContainer, ENTITY, RF> CHLView;
      mutable CHLView chLocal;

      Param& param; // parameter functions

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

    public:


      // pattern assembly flags
      enum { doPatternVolume = true };

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



      TMIXED_CHNS_NavierStokes (Param& param_, CHGFS& chGfs_, CHContainer& chV_,
                         NSGFS& nsGfs_) :
        chLocal(chGfs_,chV_),
        param(param_)
      {
      }


      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;
        //Quadrature Rule
        const int dim = EG::Entity::dimension;
        const int order = 5; //TODO: Choose more educated
        auto rule = Dune::PDELab::quadratureRule(eg.geometry(), order);

        chLocal.DoEvaluation(eg);

        for(const auto& ip : rule)
        {
          const auto S = eg.geometry().jacobianInverseTransposed(ip.position());
          auto factor = ip.weight()*eg.geometry().integrationElement(ip.position());

          CHVariables_Threephase_Micro_Extension<RF,dim,decltype(chLocal.lfs)> ch(chCache, ip.position(), chLocal.lfs, chLocal.lv, S);
          NSMicroVariables<RF,dim,LFSV> ns(vCache, ip.position(), lfsv, x.base(), S);
          RF pf3 = 1-ch.pf - ch.pf2;
          RF pf_f=ch.pf+ch.pf2+2*param.delta*pf3;

          for (size_t i=0; i<ns.vxspace.size(); i++)
          {
            r.accumulate(ns.vxspace, i, factor*(
              param.vDis.eval(pf_f)*ns.vx*ns.vbasis.phi[i]
              + param.phys.eval_mu(pf_f) * (ns.gradvx * ns.vbasis.gradphi[i][0])
              - pf_f * ns.vbasis.phi[i]
            ));
          } // for i
        } // for ip
//        std::cout << "bye_alpha" << std::endl;
      } //alpha volume

    };