Newer
Older
#ifndef DUNE_PDELAB_SYSTEMTESTING_HELMHOLTZ_PARAMETERS_SPHERICALWAVE_HH
#define DUNE_PDELAB_SYSTEMTESTING_HELMHOLTZ_PARAMETERS_SPHERICALWAVE_HH
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
/** \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;
};