dune-pdelab  2.0.0
poisson.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 
4 #ifndef DUNE_PDELAB_POISSON_HH
5 #define DUNE_PDELAB_POISSON_HH
6 
7 #include<vector>
8 
9 #include<dune/common/exceptions.hh>
10 #include<dune/common/fvector.hh>
11 #include<dune/common/static_assert.hh>
12 
13 #include<dune/geometry/type.hh>
14 #include<dune/geometry/quadraturerules.hh>
15 
16 #include <dune/localfunctions/common/interfaceswitch.hh>
17 
19 
20 #include"defaultimp.hh"
21 #include"idefault.hh"
22 #include"pattern.hh"
23 #include"flags.hh"
24 
25 namespace Dune {
26  namespace PDELab {
30 
43  template<typename F, typename B, typename J>
44  class Poisson : public NumericalJacobianApplyVolume<Poisson<F,B,J> >,
45  public FullVolumePattern,
47  {
48  public:
49  // pattern assembly flags
50  enum { doPatternVolume = true };
51 
52  // residual assembly flags
53  enum { doAlphaVolume = true };
54  enum { doLambdaVolume = true };
55  enum { doLambdaBoundary = true };
56 
61  Poisson (const F& f_, const B& bctype_, const J& j_, unsigned int quadOrder)
62  : f(f_), bctype(bctype_), j(j_),
63  laplace_(quadOrder),
64  quadOrder_(quadOrder)
65  {}
66 
67  // volume integral depending on test and ansatz functions
68  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
69  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
70  {
71  laplace_.alpha_volume(eg, lfsu, x, lfsv, r);
72  }
73 
84  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
85  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, M & matrix) const
86  {
87  laplace_.jacobian_volume(eg, lfsu, x, lfsv, matrix);
88  }
89  // volume integral depending only on test functions
90  template<typename EG, typename LFSV, typename R>
91  void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
92  {
93  // domain and range field type
94  typedef FiniteElementInterfaceSwitch<
95  typename LFSV::Traits::FiniteElementType
96  > FESwitch;
97  typedef BasisInterfaceSwitch<
98  typename FESwitch::Basis
99  > BasisSwitch;
100  typedef typename BasisSwitch::DomainField DF;
101  typedef typename BasisSwitch::RangeField RF;
102  typedef typename BasisSwitch::Range Range;
103 
104  // dimensions
105  static const int dimLocal = EG::Geometry::mydimension;
106 
107  // select quadrature rule
108  Dune::GeometryType gt = eg.geometry().type();
109  const Dune::QuadratureRule<DF,dimLocal>& rule =
110  Dune::QuadratureRules<DF,dimLocal>::rule(gt,quadOrder_);
111 
112  // loop over quadrature points
113  for (typename Dune::QuadratureRule<DF,dimLocal>::const_iterator it =
114  rule.begin(); it!=rule.end(); ++it)
115  {
116  // evaluate shape functions
117  std::vector<Range> phi(lfsv.size());
118  FESwitch::basis(lfsv.finiteElement()).
119  evaluateFunction(it->position(),phi);
120 
121  // evaluate right hand side parameter function
122  typename F::Traits::RangeType y(0.0);
123  f.evaluate(eg.entity(),it->position(),y);
124 
125  // integrate f
126  RF factor = - r.weight() * it->weight() * eg.geometry().integrationElement(it->position());
127  for (size_t i=0; i<lfsv.size(); i++)
128  r.rawAccumulate(lfsv,i,y*phi[i]*factor);
129  }
130  }
131 
132  // boundary integral independen of ansatz functions
133  template<typename IG, typename LFSV, typename R>
134  void lambda_boundary (const IG& ig, const LFSV& lfsv, R& r) const
135  {
136  // domain and range field type
137  typedef FiniteElementInterfaceSwitch<
138  typename LFSV::Traits::FiniteElementType
139  > FESwitch;
140  typedef BasisInterfaceSwitch<
141  typename FESwitch::Basis
142  > BasisSwitch;
143  typedef typename BasisSwitch::DomainField DF;
144  typedef typename BasisSwitch::DomainLocal DomainLocal;
145  typedef typename BasisSwitch::RangeField RF;
146  typedef typename BasisSwitch::Range Range;
147 
148  // dimensions
149  static const int dimLocal = IG::Geometry::mydimension;
150 
151  // select quadrature rule
152  Dune::GeometryType gtface = ig.geometryInInside().type();
153  const Dune::QuadratureRule<DF,dimLocal>& rule =
154  Dune::QuadratureRules<DF,dimLocal>::rule(gtface,quadOrder_);
155 
156  // loop over quadrature points and integrate normal flux
157  for (typename Dune::QuadratureRule<DF,dimLocal>::const_iterator it =
158  rule.begin(); it!=rule.end(); ++it)
159  {
160  // evaluate boundary condition type
161  // skip rest if we are on Dirichlet boundary
162  if( !bctype.isNeumann( ig,it->position() ) )
163  continue;
164 
165  // position of quadrature point in local coordinates of element
166  const DomainLocal& local =
167  ig.geometryInInside().global(it->position());
168 
169  // evaluate test shape functions
170  std::vector<Range> phi(lfsv.size());
171  FESwitch::basis(lfsv.finiteElement()).evaluateFunction(local,phi);
172 
173  // evaluate flux boundary condition
174  typename J::Traits::RangeType y(0.0);
175  j.evaluate(*(ig.inside()),local,y);
176 
177  // integrate J
178  RF factor = r.weight() * it->weight()*ig.geometry().integrationElement(it->position());
179  for (size_t i=0; i<lfsv.size(); i++)
180  r.rawAccumulate(lfsv,i,y*phi[i]*factor);
181  }
182  }
183 
184  protected:
185  const F& f;
186  const B& bctype;
187  const J& j;
188 
189  // Laplace assembler to handle the matrix assembly
191 
192  // Quadrature rule order
193  unsigned int quadOrder_;
194  };
195 
198 
212  template<typename Time, typename F, typename B, typename J>
214  : public Poisson<F,B,J>,
216  {
217  typedef Poisson<F,B,J> Base;
219 
220  protected:
221  // need non-const references for setTime()
222  F& f;
223  B& bctype;
224  J& j;
225 
226  public:
228  InstationaryPoisson(F& f_, B& bctype_, J& j_, unsigned int quadOrder)
229  : Base(f_, bctype_, j_, quadOrder)
230  , f(f_)
231  , bctype(bctype_)
232  , j(j_)
233  {}
234 
236  void setTime (Time t) {
237  f.setTime(t);
238  bctype.setTime(t);
239  j.setTime(t);
241  }
242  };
243 
245  } // namespace PDELab
246 } // namespace Dune
247 
248 #endif
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &matrix) const
Compute the Laplace stiffness matrix for the element given in 'eg'.
Definition: laplace.hh:126
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: poisson.hh:69
Default flags for all local operators.
Definition: flags.hh:18
const J & j
Definition: poisson.hh:187
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Compute Laplace matrix times a given vector for one element.
Definition: laplace.hh:67
const B & bctype
Definition: poisson.hh:186
unsigned int quadOrder_
Definition: poisson.hh:193
a local operator for solving the Poisson equation in instationary problems
Definition: poisson.hh:213
J & j
Definition: poisson.hh:224
const IG & ig
Definition: common/constraints.hh:146
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: poisson.hh:134
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &matrix) const
Compute the Laplace stiffness matrix for the element given in 'eg'.
Definition: poisson.hh:85
F & f
Definition: poisson.hh:222
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Poisson(const F &f_, const B &bctype_, const J &j_, unsigned int quadOrder)
Constructor.
Definition: poisson.hh:61
const F & f
Definition: poisson.hh:185
Definition: poisson.hh:44
sparsity pattern generator
Definition: pattern.hh:14
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
Definition: poisson.hh:54
const EG & eg
Definition: common/constraints.hh:277
InstationaryPoisson(F &f_, B &bctype_, J &j_, unsigned int quadOrder)
construct InstationaryPoisson
Definition: poisson.hh:228
Definition: poisson.hh:53
Laplace laplace_
Definition: poisson.hh:190
void setTime(Time t)
set the time for subsequent evaluation on the parameter functions
Definition: poisson.hh:236
B & bctype
Definition: poisson.hh:223
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: poisson.hh:91
Definition: laplace.hh:37
void setTime(Timet_)
set time for subsequent evaluation
Definition: idefault.hh:104