dgnavierstokesparameter.hh

// -*- tab-width: 2; indent-tabs-mode: nil -*-
// vi: set et ts=2 sw=2 sts=2:
 
#ifndef DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKESPARAMETER_HH
#define DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKESPARAMETER_HH
 
#include <dune/common/parametertreeparser.hh>
#include <dune/pdelab/common/geometrywrapper.hh>
#include <dune/pdelab/localoperator/dginteriorpenaltyparameter.hh>
#include <dune/pdelab/localoperator/stokesparameter.hh>
 
namespace Dune {
  namespace PDELab {
 
    template<typename GV, typename RF, typename F, typename B, typename V, typename J,
             bool navier = false, bool full_tensor = false, typename IP = DefaultInteriorPenalty<typename V::Traits::RangeFieldType> >
    class DGNavierStokesParameters :
      public NavierStokesDefaultParameters<GV,RF,F,B,V,J,navier,full_tensor>
    {
 
      typedef NavierStokesDefaultParameters<GV,RF,F,B,V,J,navier,full_tensor> Base;
 
      void initFromString(const std::string & method)
      {
        std::string s = method;
        std::transform(s.begin(), s.end(), s.begin(), tolower);
 
        // nipg (epsilon=1) 2d p1 -> Klaus sagt sollte auch sigma 1 klappen
        if (s.find("nipg") != std::string::npos)
          {
            _epsilon = 1;
            return;
          }
        // sipg (epsilon=-1) 2d p1 -> Klaus sagt sigma=3.9irgendwas
        if (s.find("sipg") != std::string::npos)
          {
            _epsilon = -1;
            return;
          }
        // obb sigma = 0, epsilon =
        if (s == "obb")
          {
            _epsilon = 1;
            return;
          }
        // extract parameters
        {
          double sigma, beta;
          if (3 == sscanf(s.c_str(), "%d %lg %lg", &_epsilon, &sigma, &beta))
            return;
        }
        DUNE_THROW(Dune::Exception, "Unknown DG type " << method);
      }
 
    public :
 
      typedef typename Base::Traits Traits;
 
      DGNavierStokesParameters(const Dune::ParameterTree& configuration,
                               F& f, B& b, V& v, J& j)
        : Base(configuration,f,b,v,j)
        , _ip(configuration.sub("dg"))
        , _epsilon(configuration.sub("dg").get<int>("epsilon"))
      {}
 
 
      typename Traits::RangeField
      incompressibilityScaling ( typename Traits::RangeField  dt ) const
      {
        typename Traits::RangeField y(1.0 / dt);
        return y;
      }
 
      template<typename GEO, typename IGEO, typename OGEO>
      typename Traits::RangeField
      getFaceIP(const GEO& geo, const IGEO& igeo, const OGEO& ogeo)
      {
        return _ip.getFaceIP(geo,igeo,ogeo);
      }
 
      template<typename GEO, typename IGEO>
      typename Traits::RangeField
      getFaceIP(const GEO& geo, const IGEO& igeo)
      {
        return _ip.getFaceIP(geo,igeo);
      }
 
      int
      epsilonIPSymmetryFactor()
      {
        return _epsilon;
      }
 
    private:
 
      IP _ip;       // Interior penalty
      int _epsilon; // IP symmetry factor
    }; // end class DGNavierStokesParameters
 
 
    namespace NavierStokesDGImp{
      template< typename PRM, typename Dummy = void>
      struct VariableBoundarySlipSwitch{
 
        template<typename IntersectionGeometry>
        static typename PRM::Traits::RangeField
        boundarySlip
        (const PRM & ,
         const IntersectionGeometry& ,
         const typename PRM::Traits::IntersectionDomain& )
        {
          return 1.0;
        }
 
      };
 
      template< typename PRM>
      struct VariableBoundarySlipSwitch
      <PRM,typename std::enable_if<PRM::enable_variable_slip>::type>
      {
        template<typename IntersectionGeometry>
        static typename PRM::Traits::RangeField
        boundarySlip
        (const PRM & prm,
         const IntersectionGeometry& ig,
         const typename PRM::Traits::IntersectionDomain& x)
        {
          return prm.boundarySlip(ig,x);
        }
 
      };
    }
 
  } // end namespace PDELab
} // end namespace Dune
#endif // DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKESPARAMETER_HH