linearacousticsparameter.hh

// -*- tab-width: 8; indent-tabs-mode: nil -*-
#ifndef DUNE_PDELAB_LOCALOPERATOR_LINEARACOUSTICSPARAMETER_HH
#define DUNE_PDELAB_LOCALOPERATOR_LINEARACOUSTICSPARAMETER_HH
 
#include<vector>
 
#include<dune/common/exceptions.hh>
#include<dune/common/fvector.hh>
#include<dune/geometry/type.hh>
#include<dune/geometry/referenceelements.hh>
#include<dune/geometry/quadraturerules.hh>
#include<dune/pdelab/common/geometrywrapper.hh>
#include<dune/pdelab/common/function.hh>
#include<dune/pdelab/common/functionutilities.hh>
 
namespace Dune {
  namespace PDELab {
 
    template<typename GV, typename RF>
    struct LinearAcousticsParameterTraits
    {
      typedef GV GridViewType;
 
      enum {
        dimDomain = GV::dimension
      };
 
      typedef typename GV::Grid::ctype DomainFieldType;
 
      typedef Dune::FieldVector<DomainFieldType,dimDomain> DomainType;
 
      typedef Dune::FieldVector<DomainFieldType,dimDomain-1> IntersectionDomainType;
 
      typedef RF RangeFieldType;
 
      typedef Dune::FieldVector<RF,GV::dimension> RangeType;
 
      typedef Dune::FieldVector<RF,GV::dimension+1> StateType;
 
      typedef typename GV::Traits::template Codim<0>::Entity ElementType;
      typedef typename GV::Intersection IntersectionType;
    };
 
    template<typename T>
    class LinearAcousticsInitialValueAdapter
      : public Dune::PDELab::GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename T::Traits::GridViewType,
                                                                               typename T::Traits::RangeFieldType,
                                                                               T::Traits::dimDomain+1,Dune::FieldVector<typename T::Traits::RangeFieldType,T::Traits::dimDomain+1> >
                                              ,LinearAcousticsInitialValueAdapter<T> >
    {
    public:
      typedef Dune::PDELab::GridFunctionTraits<typename T::Traits::GridViewType,
                                               typename T::Traits::RangeFieldType,
                                               T::Traits::dimDomain+1,Dune::FieldVector<typename T::Traits::RangeFieldType,T::Traits::dimDomain+1> > Traits;
 
      LinearAcousticsInitialValueAdapter (const typename Traits::GridViewType& g_, const T& t_) : g(g_), t(t_) {}
 
      inline void evaluate (const typename Traits::ElementType& e,
                            const typename Traits::DomainType& x,
                            typename Traits::RangeType& y) const
      {
        y = t.u0(e,x);
      }
 
      inline const typename Traits::GridViewType& getGridView () const
      {
        return g;
      }
 
    private:
      typename Traits::GridViewType g;
      const T& t;
    };
 
  template<typename GV, typename RF>
  class LinearAcousticsModelProblem
  {
  public:
    typedef LinearAcousticsParameterTraits<GV,RF> Traits;
 
    LinearAcousticsModelProblem ()
      : pi(3.141592653589793238462643), time(0.0)
    {
    }
 
    typename Traits::RangeFieldType
    c (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
    {
      return 340.0;
    }
 
    typename Traits::StateType
    g (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x, const typename Traits::StateType& s) const
    {
      typename Traits::DomainType xglobal = is.geometry().global(x);
      if (xglobal[0]<1e-6)
        {
          typename Traits::StateType u(0.0);
          u[0] = s[0];
          u[1] = 1.224*(1+0.5*sin(2*pi*1500.0*time));
          return u;
        }
      if (xglobal[0]>1.0-1e-6)
        {
          typename Traits::StateType u(0.0);
          return u;
        }
      typename Traits::StateType u(0.0);
      u[0] = s[0];
      u[2] = 0.0;
      return u;
    }
 
    typename Traits::StateType
    q (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
    {
      typename Traits::StateType rhs(0.0);
      return rhs;
    }
 
    typename Traits::StateType
    u0 (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
    {
      typename Traits::StateType u(0.0);
      return u;
    }
 
    void setTime (RF t)
    {
      time = t;
    }
 
  private:
    double pi;
    RF time;
  };
 
 
 
  }
}
 
#endif // DUNE_PDELAB_LOCALOPERATOR_LINEARACOUSTICSPARAMETER_HH