// -*- tab-width: 4; indent-tabs-mode: nil -*- /** \file \brief Plane Wave Parameter = Problem description */ #ifndef PARAMETERSPLANEWAVE_HH #define PARAMETERSPLANEWAVE_HH /** \brief Defines problem parameter analytical solution */ template<typename GV, typename RF, typename DF> class ParametersPlaneWave : public Dune::PDELab::DirichletConstraintsParameters { public: typedef Dune::PDELab::GridFunctionTraits<GV,RF,1,Dune::FieldVector<RF,1> > Traits; // typedef DiffusionParameterTraits<GV,RF> Traits; ParametersPlaneWave (double omega_) : omega(omega_) { } //! Dirichlet boundary condition type function template<typename IG> inline bool isDirichlet(const IG& intersection, const Dune::FieldVector<typename IG::ctype, IG::dimension-1>& xlocal) const { return false; } //! Neumann boundary condition type function template<typename IG> inline bool isNeumann(const IG& intersection, const Dune::FieldVector<typename IG::ctype, IG::dimension-1>& xlocal) const { return true; } //! Dirichlet boundary condition value RF g (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal ) const { return RF(0.,0.); } //! Neumann boundary condition value template <typename IG> RF j (const IG& ig, const Dune::FieldVector<typename IG::ctype, IG::dimension-1>& xlocal, RF u) const { const int dim = Traits::GridViewType::Grid::dimension; typedef typename Traits::GridViewType::Grid::ctype ctype; Dune::FieldVector<ctype,dim> xglobal = ig.geometry().global(xlocal); RF i(0., 1.); RF x,y; x = xglobal[0]; y = xglobal[1]; // get the outer normal vector at the face center const Dune::FieldVector<RF, dim> normal = ig.centerUnitOuterNormal(); Dune::FieldVector<RF, dim> gradu; gradu[0] = i*omega*cos(theta)*std::exp( i*omega*(x*cos(theta) + y *sin(theta)) ); gradu[1] = i*omega*sin(theta)*std::exp( i*omega*(x*cos(theta) + y *sin(theta)) ); RF res = gradu * normal; res *= -1; return res; } //! source term template <typename EG> RF f (const EG& eg, const Dune::FieldVector<typename EG::Geometry::ctype, EG::Geometry::dimension>& position) const { RF res(0.,0.); return res; } //! exact solution / analytic solution RF ua (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal ) const { const int dim = Traits::GridViewType::Grid::dimension; typedef typename Traits::GridViewType::Grid::ctype ctype; Dune::FieldVector<ctype,dim> xglobal = e.geometry().global(xlocal); RF i(0., 1.); RF x,y; x = xglobal[0]; y = xglobal[1]; return std::exp( i*omega*(x*cos(theta) + y *sin(theta)) ); } const static int dim = Traits::GridViewType::Grid::dimension; //! exact grad of solution / analytic grad of solution Dune::FieldVector<RF,dim> ugrada (const typename Traits::ElementType& e, const typename Traits::DomainType& xlocal ) const { typedef typename Traits::GridViewType::Grid::ctype ctype; Dune::FieldVector<ctype,dim> xglobal = e.geometry().global(xlocal); RF i(0., 1.); RF x,y; x = xglobal[0]; y = xglobal[1]; Dune::FieldVector<RF,dim> res; res[0] = i*omega*cos(theta)*std::exp( i*omega*(x*cos(theta) + y *sin(theta)) ); res[1] = i*omega*sin(theta)*std::exp( i*omega*(x*cos(theta) + y *sin(theta)) ); return res; } double theta = M_PI/2.; double omega; }; #endif // PARAMETERSPLANEWAVE_HH