stokesparameter.hh

#ifndef DUNE_PDELAB_LOCALOPERATOR_STOKESPARAMETER_HH
#define DUNE_PDELAB_LOCALOPERATOR_STOKESPARAMETER_HH
 
#include <dune/common/parametertree.hh>
#include <dune/pdelab/common/function.hh>
#include <dune/pdelab/constraints/common/constraintsparameters.hh>
 
namespace Dune {
  namespace PDELab {
 
    struct StokesBoundaryCondition {
      enum Type {
        DoNothing = 0,
        VelocityDirichlet = 1,
        StressNeumann = 2,
        SlipVelocity = 3
      };
    };
 
    template<typename GV, typename RF>
    struct NavierStokesParameterTraits
    {
      typedef GV GridView;
 
      enum {
        dimDomain = GV::dimension
      };
 
      typedef typename GV::Grid::ctype DomainField;
 
      typedef Dune::FieldVector<DomainField,dimDomain> Domain;
 
      typedef Dune::FieldVector<DomainField,dimDomain-1> IntersectionDomain;
 
      typedef RF RangeField;
 
      typedef Dune::FieldVector<RF,GV::dimensionworld> VelocityRange;
 
      typedef Dune::FieldVector<RF,1> PressureRange;
 
      typedef StokesBoundaryCondition BoundaryCondition;
 
      typedef typename GV::Traits::template Codim<0>::Entity Element;
      typedef typename GV::Intersection Intersection;
    };
 
    namespace {
 
      template<typename GF, typename Entity, typename Domain>
      typename GF::Traits::RangeType
      evaluateVelocityGridFunction(const GF& gf,
                                   const Entity& e,
                                   const Domain& x)
      {
        static_assert(int(GF::Traits::dimRange) == int(Domain::dimension),"dimension of function range does not match grid dimension");
        typename GF::Traits::RangeType y;
        gf.evaluate(e,x,y);
        return y;
      }
 
      template<typename GF, typename Entity, typename Domain>
      FieldVector<typename GF::template Child<0>::Type::Traits::RangeFieldType,TypeTree::StaticDegree<GF>::value>
        evaluateVelocityGridFunction(const GF& gf,
                                     const Entity& e,
                                     const Domain& x)
      {
        static_assert(Domain::dimension == TypeTree::StaticDegree<GF>::value,"dimension of function range does not match grid dimension");
        FieldVector<typename GF::template Child<0>::Type::Traits::RangeFieldType,TypeTree::StaticDegree<GF>::value> y;
        typename GF::template Child<0>::Type::Traits::RangeType cy;
        for (int d = 0; d < Domain::dimension; ++d)
          {
            gf.child(d).evaluate(e,x,cy);
            y[d] = cy;
          }
        return y;
      }
 
    }
 
    template <typename GV, typename RF, typename F, typename B, typename V, typename J, bool navier = false, bool tensor = false>
    class NavierStokesDefaultParameters
    {
    public:
 
      static const bool assemble_navier = navier;
      static const bool assemble_full_tensor = tensor;
 
      typedef NavierStokesParameterTraits<GV,RF> Traits;
 
      NavierStokesDefaultParameters(const Dune::ParameterTree& config,
                                    F& f,
                                    B& b,
                                    V& v,
                                    J& j)
        : _rho(config.get<RF>("rho"))
        , _mu(config.get<RF>("mu"))
        , _f(f)
        , _b(b)
        , _v(v)
        , _j(j)
      {}
 
      NavierStokesDefaultParameters(const RF& mu,
                                    const RF& rho,
                                    F& f,
                                    B& b,
                                    V& v,
                                    J& j)
        : _rho(rho)
        , _mu(mu)
        , _f(f)
        , _b(b)
        , _v(v)
        , _j(j)
      {}
 
 
      template<typename EG>
      typename Traits::VelocityRange
      f(const EG& e, const typename Traits::Domain& x) const
      {
        typename F::Traits::RangeType fvalue;
        return evaluateVelocityGridFunction(_f,e,x);
      }
 
      template<typename IG>
      typename Traits::BoundaryCondition::Type
      bctype(const IG& is,
             const typename Traits::IntersectionDomain& x) const
      {
        typename B::Traits::BoundaryCondition::Type y;
        _b.evaluate(is,x,y);
        return y;
      }
 
      template<typename EG>
      typename Traits::RangeField
      mu(const EG& e, const typename Traits::Domain& x) const
      {
        return _mu;
      }
 
      template<typename IG>
      typename Traits::RangeField
      mu(const IG& ig, const typename Traits::IntersectionDomain& x) const
      {
        return _mu;
      }
 
      template<typename EG>
      typename Traits::RangeField
      rho(const EG& eg, const typename Traits::Domain& x) const
      {
        return _rho;
      }
 
      template<typename IG>
      typename Traits::RangeField
      rho(const IG& ig, const typename Traits::IntersectionDomain& x) const
      {
        return _rho;
      }
 
      template<typename EG>
      typename Traits::VelocityRange
      g(const EG& e, const typename Traits::Domain& x) const
      {
        typename V::Traits::RangeType y;
        _v.evaluate(e,x,y);
        return y;
      }
 
      template<typename EG>
      typename Traits::RangeField
      g2(const EG& e, const typename Traits::Domain& x) const
      {
        return 0;
      }
 
#ifdef DOXYGEN
 
      template<typename IG>
      typename Traits::VelocityRange>
      j(const IG& ig,
        const typename Traits::IntersectionDomain& x,
        const typename Traits::Domain& normal) const;
 
#else // DOXYGEN
 
      template<typename IG>
      typename std::enable_if<
        J::Traits::dimRange == 1 &&
        (GV::dimension > 1) &&
        AlwaysTrue<IG>::value, // required to force lazy evaluation
        typename Traits::VelocityRange
          >::type
          j(const IG& ig,
            const typename Traits::IntersectionDomain& x,
            typename Traits::Domain normal) const
      {
        typename J::Traits::RangeType r;
        auto e = ig.inside();
        _j.evaluate(e,ig.geometryInInside().global(x),r);
        normal *= r;
        return normal;
      }
 
      template<typename IG>
      typename std::enable_if<
        J::Traits::dimRange == GV::dimension &&
        AlwaysTrue<IG>::value, // required to force lazy evaluation
        typename Traits::VelocityRange
        >::type
        j(const IG& ig,
          const typename Traits::IntersectionDomain& x,
          const typename Traits::Domain& normal) const
      {
        auto e = ig.inside();
        typename J::Traits::RangeType y;
        _j.evaluate(e,ig.geometryInInside().global(x),y);
        return y;
      }
 
#endif // DOXYGEN
 
      void setTime(RF time)
      {
        _f.setTime(time);
        _b.setTime(time);
        _v.setTime(time);
        _j.setTime(time);
      }
 
    private:
      const RF _rho;
      const RF _mu;
      F& _f;
      B& _b;
      V& _v;
      J& _j;
    };
 
 
    template<typename PRM>
    class StokesVelocityDirichletConstraints
      : public Dune::PDELab::DirichletConstraintsParameters
    {
    private:
      const PRM & prm_;
 
    public:
 
      StokesVelocityDirichletConstraints (const PRM & _prm)
        : prm_(_prm) { }
 
      template<typename I>
      bool isDirichlet(const I & intersection,
                       const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
      {
        StokesBoundaryCondition::Type bctype = prm_.bctype(intersection,coord);
        return (bctype == StokesBoundaryCondition::VelocityDirichlet);
      }
 
    };
 
    template<typename PRM>
    class StokesPressureDirichletConstraints
      : public Dune::PDELab::DirichletConstraintsParameters
    {
    private:
      const PRM & prm_;
 
    public:
 
      StokesPressureDirichletConstraints (const PRM & _prm)
        : prm_(_prm) { }
 
      template<typename I>
      bool isDirichlet(const I & intersection,
                       const Dune::FieldVector<typename I::ctype, I::mydimension> & coord) const
      { return false; }
    };
 
 
 
#ifndef DOXYGEN
 
    template<typename PRM, int rangeDim>
    class NavierStokesFunctionAdapterBase
      : public Dune::PDELab::GridFunctionBase<
      Dune::PDELab::GridFunctionTraits<
        typename PRM::Traits::GridView,
        typename PRM::Traits::RangeField,
        rangeDim,
        Dune::FieldVector<typename PRM::Traits::RangeField,rangeDim>
        >,
      NavierStokesFunctionAdapterBase<PRM,rangeDim>
      >
    {
    public:
      typedef Dune::PDELab::GridFunctionTraits<
      typename PRM::Traits::GridView,
      typename PRM::Traits::RangeField,
      rangeDim,
      Dune::FieldVector<typename PRM::Traits::RangeField,rangeDim>
      > Traits;
 
      NavierStokesFunctionAdapterBase(PRM& prm)
        : _prm(prm)
      {}
 
      void setTime(const double time)
      {
        _prm.setTime(time);
      }
 
      const PRM& parameters() const
      {
        return _prm;
      }
 
      const typename Traits::GridViewType& getGridView () const
      {
        return _prm.gridView();
      }
 
    private:
      PRM& _prm;
    };
 
 
#endif // DOXYGEN
 
    template<typename PRM>
    class NavierStokesVelocityFunctionAdapter
      : public NavierStokesFunctionAdapterBase<PRM,PRM::Traits::dimDomain>
    {
 
      typedef NavierStokesFunctionAdapterBase<PRM,PRM::Traits::dimDomain> Base;
 
      using Base::parameters;
 
    public:
 
      typedef typename Base::Traits Traits;
 
      NavierStokesVelocityFunctionAdapter(PRM& prm)
        : Base(prm)
      {}
 
      void evaluate (const typename Traits::ElementType& e,
                     const typename Traits::DomainType& x,
                     typename Traits::RangeType& y) const
      {
        y = parameters().g(e,x);
      }
    };
 
 
 
#if 0
    template < typename PRM >
    class NavierStokesDirichletFunctionAdapterFactory
    {
    public:
      typedef Dune::PDELab::CompositeGridFunction<
      NavierStokesDirichletFunctionAdapter<PRM>,
      NavierStokesPressureDirichletFunctionAdapter<PRM> >
      BoundaryDirichletFunction;
 
      NavierStokesDirichletFunctionAdapterFactory(PRM & prm)
        : v(prm), p(prm), df(v,p)
      {}
 
      BoundaryDirichletFunction & dirichletFunction()
      {
        return df;
      }
 
    private:
      NavierStokesVelocityDirichletFunctionAdapter<PRM> v;
      NavierStokesPressureDirichletFunctionAdapter<PRM> p;
      BoundaryDirichletFunction df;
    };
#endif
 
 
  } // end namespace PDELab
} // end namespace Dune
 
#endif // DUNE_PDELAB_LOCALOPERATOR_STOKESPARAMETER_HH