// -*- tab-width: 4; indent-tabs-mode: nil -*-

#ifndef DUNE_PDELAB_SYSTEMTESTING_HELMHOLTZ_PARAMETERS_SPHERICALWAVE_HH
#define DUNE_PDELAB_SYSTEMTESTING_HELMHOLTZ_PARAMETERS_SPHERICALWAVE_HH


/** \file

    \brief Spherical Wave Parameter = Problem description
*/
/** \brief Defines problem parameter: Dirichlet boundary conditions AND its
 *         extension to the interior, sinks and sources and the reaction rate.
 */
template<typename GV, typename RF, typename DF>
class ParametersSphericalWave  : public Dune::PDELab::DirichletConstraintsParameters
{
public:

  typedef Dune::PDELab::GridFunctionTraits<GV,RF,1,Dune::FieldVector<RF,1> > Traits;
  // typedef DiffusionParameterTraits<GV,RF> Traits;

  ParametersSphericalWave (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
  {
    RF i(0., 1.);
    return i*omega*u;
  }



    //! source term
    template <typename EG>
    RF f (const EG& eg,  const Dune::FieldVector<typename EG::Geometry::ctype, EG::Geometry::dimension>& position) const
    {


        Dune::FieldVector<typename EG::Geometry::ctype, EG::Geometry::dimension> xg = eg.geometry().global( position );
        Dune::FieldVector<DF,EG::Geometry::dimension> midpoint;

        for(int i=0; i<EG::Geometry::dimension; ++i) {
            midpoint[i] = 0.5;
        }

        xg -= midpoint;

        RF res(0.,0.);

        double r=0.2;

        double tmp=0.;
        for(int i=0; i<EG::Geometry::dimension; ++i) {
            tmp += xg[i]*xg[i];
        }

        if( std::sqrt( tmp ) < r )
            res = RF(25,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);

        std::cout<<"Error: Analytical Solution has not been implemented"<<std::endl;

        return RF(0.,0.);

    }

    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);

        std::cout<<"Error: Analytical Solution has not been implemented"<<std::endl;
        Dune::FieldVector<RF,dim> res;
        res[0] = 0.;
        res[1] = 0.;
        return res;

    }


  double omega;

};

#endif