dune-pdelab  2.0.0
convectiondiffusiondg.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_CONVECTIONDIFFUSIONDG_HH
3 #define DUNE_PDELAB_CONVECTIONDIFFUSIONDG_HH
4 
5 #include<dune/common/exceptions.hh>
6 #include<dune/common/fvector.hh>
7 #include<dune/common/static_assert.hh>
8 #include<dune/geometry/referenceelements.hh>
15 
17 
18 #ifndef USECACHE
19 #define USECACHE 1
20 //#define USECACHE 0
21 #endif
22 
30 namespace Dune {
31  namespace PDELab {
32 
34  {
35  enum Type { NIPG, SIPG };
36  };
37 
39  {
41  };
42 
57  template<typename T, typename FiniteElementMap>
59  : public Dune::PDELab::NumericalJacobianApplyVolume<ConvectionDiffusionDG<T,FiniteElementMap> >,
60  public Dune::PDELab::NumericalJacobianApplySkeleton<ConvectionDiffusionDG<T,FiniteElementMap> >,
61  public Dune::PDELab::NumericalJacobianApplyBoundary<ConvectionDiffusionDG<T,FiniteElementMap> >,
65  public Dune::PDELab::InstationaryLocalOperatorDefaultMethods<typename T::Traits::RangeFieldType>
66  {
67  enum { dim = T::Traits::GridViewType::dimension };
68 
69  typedef typename T::Traits::RangeFieldType Real;
71 
72  public:
73  // pattern assembly flags
74  enum { doPatternVolume = true };
75  enum { doPatternSkeleton = true };
76 
77  // residual assembly flags
78  enum { doAlphaVolume = true };
79  enum { doAlphaSkeleton = true };
80  enum { doAlphaBoundary = true };
81  enum { doLambdaVolume = true };
82 
87  Real alpha_=0.0,
88  int intorderadd_=0
89  )
90  : Dune::PDELab::NumericalJacobianApplyVolume<ConvectionDiffusionDG<T,FiniteElementMap> >(1.0e-7),
91  Dune::PDELab::NumericalJacobianApplySkeleton<ConvectionDiffusionDG<T,FiniteElementMap> >(1.0e-7),
92  Dune::PDELab::NumericalJacobianApplyBoundary<ConvectionDiffusionDG<T,FiniteElementMap> >(1.0e-7),
93  param(param_), method(method_), weights(weights_),
94  alpha(alpha_), intorderadd(intorderadd_), quadrature_factor(2),
95  cache(20)
96  {
97  theta = 1.0;
98  if (method==ConvectionDiffusionDGMethod::SIPG) theta = -1.0;
99  }
100 
101  // volume integral depending on test and ansatz functions
102  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
103  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
104  {
105  // domain and range field type
106  typedef typename LFSU::Traits::FiniteElementType::
107  Traits::LocalBasisType::Traits::DomainFieldType DF;
108  typedef typename LFSU::Traits::FiniteElementType::
109  Traits::LocalBasisType::Traits::RangeFieldType RF;
110  typedef typename LFSU::Traits::FiniteElementType::
111  Traits::LocalBasisType::Traits::JacobianType JacobianType;
112  typedef typename LFSU::Traits::FiniteElementType::
113  Traits::LocalBasisType::Traits::RangeType RangeType;
114  typedef typename LFSU::Traits::SizeType size_type;
115 
116  // dimensions
117  const int dim = EG::Geometry::dimension;
118  const int order = std::max(lfsu.finiteElement().localBasis().order(),
119  lfsv.finiteElement().localBasis().order());
120  const int intorder = intorderadd + quadrature_factor * order;
121 
122  // select quadrature rule
123  Dune::GeometryType gt = eg.geometry().type();
124  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
125 
126  // evaluate diffusion tensor at cell center, assume it is constant over elements
127  typename T::Traits::PermTensorType A;
128  Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
129  A = param.A(eg.entity(),localcenter);
130 
131  // transformation
132  typename EG::Geometry::JacobianInverseTransposed jac;
133 
134  // loop over quadrature points
135  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
136  {
137  // evaluate basis functions
138 #if USECACHE==0
139  std::vector<RangeType> phi(lfsu.size());
140  lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
141  std::vector<RangeType> psi(lfsv.size());
142  lfsv.finiteElement().localBasis().evaluateFunction(it->position(),psi);
143 #else
144  const std::vector<RangeType>& phi = cache[order].evaluateFunction(it->position(),lfsu.finiteElement().localBasis());
145  const std::vector<RangeType>& psi = cache[order].evaluateFunction(it->position(),lfsv.finiteElement().localBasis());
146 #endif
147 
148  // evaluate u
149  RF u=0.0;
150  for (size_type i=0; i<lfsu.size(); i++)
151  u += x(lfsu,i)*phi[i];
152 
153  // evaluate gradient of basis functions (we assume Galerkin method lfsu=lfsv)
154 #if USECACHE==0
155  std::vector<JacobianType> js(lfsu.size());
156  lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
157  std::vector<JacobianType> js_v(lfsv.size());
158  lfsv.finiteElement().localBasis().evaluateJacobian(it->position(),js_v);
159 #else
160  const std::vector<JacobianType>& js = cache[order].evaluateJacobian(it->position(),lfsu.finiteElement().localBasis());
161  const std::vector<JacobianType>& js_v = cache[order].evaluateJacobian(it->position(),lfsv.finiteElement().localBasis());
162 #endif
163 
164  // transform gradients of shape functions to real element
165  jac = eg.geometry().jacobianInverseTransposed(it->position());
166  std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
167  for (size_type i=0; i<lfsu.size(); i++)
168  jac.mv(js[i][0],gradphi[i]);
169 
170  std::vector<Dune::FieldVector<RF,dim> > gradpsi(lfsv.size());
171  for (size_type i=0; i<lfsv.size(); i++)
172  jac.mv(js_v[i][0],gradpsi[i]);
173 
174  // compute gradient of u
175  Dune::FieldVector<RF,dim> gradu(0.0);
176  for (size_type i=0; i<lfsu.size(); i++)
177  gradu.axpy(x(lfsu,i),gradphi[i]);
178 
179  // compute A * gradient of u
180  Dune::FieldVector<RF,dim> Agradu(0.0);
181  A.umv(gradu,Agradu);
182 
183  // evaluate velocity field
184  typename T::Traits::RangeType b = param.b(eg.entity(),it->position());
185 
186  // evaluate reaction term
187  typename T::Traits::RangeFieldType c = param.c(eg.entity(),it->position());
188 
189  // integrate (A grad u - bu)*grad phi_i + a*u*phi_i
190  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
191  for (size_type i=0; i<lfsv.size(); i++)
192  r.accumulate(lfsv,i,( Agradu*gradpsi[i] - u*(b*gradpsi[i]) + c*u*psi[i] )*factor);
193  }
194  }
195 
196  // jacobian of volume term
197  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
198  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
199  M& mat) const
200  {
201  // domain and range field type
202  typedef typename LFSU::Traits::FiniteElementType::
203  Traits::LocalBasisType::Traits::DomainFieldType DF;
204  typedef typename LFSU::Traits::FiniteElementType::
205  Traits::LocalBasisType::Traits::RangeFieldType RF;
206  typedef typename LFSU::Traits::FiniteElementType::
207  Traits::LocalBasisType::Traits::JacobianType JacobianType;
208  typedef typename LFSU::Traits::FiniteElementType::
209  Traits::LocalBasisType::Traits::RangeType RangeType;
210  typedef typename LFSU::Traits::SizeType size_type;
211 
212  // dimensions
213  const int dim = EG::Geometry::dimension;
214  const int order = std::max(lfsu.finiteElement().localBasis().order(),
215  lfsv.finiteElement().localBasis().order());
216  const int intorder = intorderadd + quadrature_factor * order;
217 
218  // select quadrature rule
219  Dune::GeometryType gt = eg.geometry().type();
220  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
221 
222  // evaluate diffusion tensor at cell center, assume it is constant over elements
223  typename T::Traits::PermTensorType A;
224  Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
225  A = param.A(eg.entity(),localcenter);
226 
227  // transformation
228  typename EG::Geometry::JacobianInverseTransposed jac;
229 
230  // loop over quadrature points
231  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
232  {
233  // evaluate basis functions
234 #if USECACHE==0
235  std::vector<RangeType> phi(lfsu.size());
236  lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
237 #else
238  const std::vector<RangeType>& phi = cache[order].evaluateFunction(it->position(),lfsu.finiteElement().localBasis());
239 #endif
240 
241  // evaluate gradient of basis functions (we assume Galerkin method lfsu=lfsv)
242 #if USECACHE==0
243  std::vector<JacobianType> js(lfsu.size());
244  lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
245 #else
246  const std::vector<JacobianType>& js = cache[order].evaluateJacobian(it->position(),lfsu.finiteElement().localBasis());
247 #endif
248 
249  // transform gradients of shape functions to real element
250  jac = eg.geometry().jacobianInverseTransposed(it->position());
251  std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
252  std::vector<Dune::FieldVector<RF,dim> > Agradphi(lfsu.size());
253  for (size_type i=0; i<lfsu.size(); i++)
254  {
255  jac.mv(js[i][0],gradphi[i]);
256  A.mv(gradphi[i],Agradphi[i]);
257  }
258 
259  // evaluate velocity field
260  typename T::Traits::RangeType b = param.b(eg.entity(),it->position());
261 
262  // evaluate reaction term
263  typename T::Traits::RangeFieldType c = param.c(eg.entity(),it->position());
264 
265  // integrate (A grad u - bu)*grad phi_i + a*u*phi_i
266  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
267  for (size_type j=0; j<lfsu.size(); j++)
268  for (size_type i=0; i<lfsu.size(); i++)
269  mat.accumulate(lfsu,i,lfsu,j,( Agradphi[j]*gradphi[i] - phi[j]*(b*gradphi[i]) + c*phi[j]*phi[i] )*factor);
270  }
271  }
272 
273  // skeleton integral depending on test and ansatz functions
274  // each face is only visited ONCE!
275  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
276  void alpha_skeleton (const IG& ig,
277  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
278  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
279  R& r_s, R& r_n) const
280  {
281  // domain and range field type
282  typedef typename LFSV::Traits::FiniteElementType::
283  Traits::LocalBasisType::Traits::DomainFieldType DF;
284  typedef typename LFSV::Traits::FiniteElementType::
285  Traits::LocalBasisType::Traits::RangeFieldType RF;
286  typedef typename LFSV::Traits::FiniteElementType::
287  Traits::LocalBasisType::Traits::RangeType RangeType;
288  typedef typename LFSU::Traits::FiniteElementType::
289  Traits::LocalBasisType::Traits::JacobianType JacobianType;
290  typedef typename LFSV::Traits::SizeType size_type;
291 
292  // dimensions
293  const int dim = IG::dimension;
294  const int order = std::max(
295  std::max(lfsu_s.finiteElement().localBasis().order(),
296  lfsu_n.finiteElement().localBasis().order()),
297  std::max(lfsv_s.finiteElement().localBasis().order(),
298  lfsv_n.finiteElement().localBasis().order())
299  );
300  const int intorder = intorderadd+quadrature_factor*order;
301 
302  // evaluate permeability tensors
303  const Dune::FieldVector<DF,dim>&
304  inside_local = Dune::ReferenceElements<DF,dim>::general(ig.inside()->type()).position(0,0);
305  const Dune::FieldVector<DF,dim>&
306  outside_local = Dune::ReferenceElements<DF,dim>::general(ig.outside()->type()).position(0,0);
307  typename T::Traits::PermTensorType A_s, A_n;
308  A_s = param.A(*(ig.inside()),inside_local);
309  A_n = param.A(*(ig.outside()),outside_local);
310 
311  // face diameter; this should be revised for anisotropic meshes?
312  DF h_s, h_n;
313  DF hmax_s = 0.;
314  DF hmax_n = 0.;
315  element_size(ig.inside()->geometry(),h_s,hmax_s);
316  element_size(ig.outside()->geometry(),h_n,hmax_n);
317  RF h_F = std::min(h_s,h_n);
318  h_F = std::min(ig.inside()->geometry().volume(),ig.outside()->geometry().volume())/ig.geometry().volume(); // Houston!
319 
320  // select quadrature rule
321  Dune::GeometryType gtface = ig.geometryInInside().type();
322  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
323 
324  // transformation
325  typename IG::Entity::Geometry::JacobianInverseTransposed jac;
326 
327  // tensor times normal
328  const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
329  Dune::FieldVector<RF,dim> An_F_s;
330  A_s.mv(n_F,An_F_s);
331  Dune::FieldVector<RF,dim> An_F_n;
332  A_n.mv(n_F,An_F_n);
333 
334  // compute weights
335  RF omega_s;
336  RF omega_n;
337  RF harmonic_average(0.0);
339  {
340  RF delta_s = (An_F_s*n_F);
341  RF delta_n = (An_F_n*n_F);
342  omega_s = delta_n/(delta_s+delta_n+1e-20);
343  omega_n = delta_s/(delta_s+delta_n+1e-20);
344  harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1e-20);
345  }
346  else
347  {
348  omega_s = omega_n = 0.5;
349  harmonic_average = 1.0;
350  }
351 
352  // get polynomial degree
353  const int order_s = lfsu_s.finiteElement().localBasis().order();
354  const int order_n = lfsu_n.finiteElement().localBasis().order();
355  int degree = std::max( order_s, order_n );
356 
357  // penalty factor
358  RF penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
359 
360  // loop over quadrature points
361  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
362  {
363  // exact normal
364  const Dune::FieldVector<DF,dim> n_F_local = ig.unitOuterNormal(it->position());
365 
366  // position of quadrature point in local coordinates of elements
367  Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(it->position());
368  Dune::FieldVector<DF,dim> iplocal_n = ig.geometryInOutside().global(it->position());
369 
370  // evaluate basis functions
371 #if USECACHE==0
372  std::vector<RangeType> phi_s(lfsu_s.size());
373  lfsu_s.finiteElement().localBasis().evaluateFunction(iplocal_s,phi_s);
374  std::vector<RangeType> phi_n(lfsu_n.size());
375  lfsu_n.finiteElement().localBasis().evaluateFunction(iplocal_n,phi_n);
376  std::vector<RangeType> psi_s(lfsv_s.size());
377  lfsv_s.finiteElement().localBasis().evaluateFunction(iplocal_s,psi_s);
378  std::vector<RangeType> psi_n(lfsv_n.size());
379  lfsv_n.finiteElement().localBasis().evaluateFunction(iplocal_n,psi_n);
380 #else
381  const std::vector<RangeType>& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
382  const std::vector<RangeType>& phi_n = cache[order_n].evaluateFunction(iplocal_n,lfsu_n.finiteElement().localBasis());
383  const std::vector<RangeType>& psi_s = cache[order_s].evaluateFunction(iplocal_s,lfsv_s.finiteElement().localBasis());
384  const std::vector<RangeType>& psi_n = cache[order_n].evaluateFunction(iplocal_n,lfsv_n.finiteElement().localBasis());
385 #endif
386 
387  // evaluate u
388  RF u_s=0.0;
389  for (size_type i=0; i<lfsu_s.size(); i++)
390  u_s += x_s(lfsu_s,i)*phi_s[i];
391  RF u_n=0.0;
392  for (size_type i=0; i<lfsu_n.size(); i++)
393  u_n += x_n(lfsu_n,i)*phi_n[i];
394 
395  // evaluate gradient of basis functions (we assume Galerkin method lfsu=lfsv)
396 #if USECACHE==0
397  std::vector<JacobianType> gradphi_s(lfsu_s.size());
398  lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
399  std::vector<JacobianType> gradphi_n(lfsu_n.size());
400  lfsu_n.finiteElement().localBasis().evaluateJacobian(iplocal_n,gradphi_n);
401  std::vector<JacobianType> gradpsi_s(lfsv_s.size());
402  lfsv_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradpsi_s);
403  std::vector<JacobianType> gradpsi_n(lfsv_n.size());
404  lfsv_n.finiteElement().localBasis().evaluateJacobian(iplocal_n,gradpsi_n);
405 #else
406  const std::vector<JacobianType>& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
407  const std::vector<JacobianType>& gradphi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsu_n.finiteElement().localBasis());
408  const std::vector<JacobianType>& gradpsi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsv_s.finiteElement().localBasis());
409  const std::vector<JacobianType>& gradpsi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsv_n.finiteElement().localBasis());
410 #endif
411 
412  // transform gradients of shape functions to real element
413  jac = ig.inside()->geometry().jacobianInverseTransposed(iplocal_s);
414  std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
415  for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
416  std::vector<Dune::FieldVector<RF,dim> > tgradpsi_s(lfsv_s.size());
417  for (size_type i=0; i<lfsv_s.size(); i++) jac.mv(gradpsi_s[i][0],tgradpsi_s[i]);
418  jac = ig.outside()->geometry().jacobianInverseTransposed(iplocal_n);
419  std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
420  for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
421  std::vector<Dune::FieldVector<RF,dim> > tgradpsi_n(lfsv_n.size());
422  for (size_type i=0; i<lfsv_n.size(); i++) jac.mv(gradpsi_n[i][0],tgradpsi_n[i]);
423 
424  // compute gradient of u
425  Dune::FieldVector<RF,dim> gradu_s(0.0);
426  for (size_type i=0; i<lfsu_s.size(); i++)
427  gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
428  Dune::FieldVector<RF,dim> gradu_n(0.0);
429  for (size_type i=0; i<lfsu_n.size(); i++)
430  gradu_n.axpy(x_n(lfsu_n,i),tgradphi_n[i]);
431 
432  // evaluate velocity field and upwinding, assume H(div) velocity field => may choose any side
433  typename T::Traits::RangeType b = param.b(*(ig.inside()),iplocal_s);
434  RF normalflux = b*n_F_local;
435  RF omegaup_s, omegaup_n;
436  if (normalflux>=0.0)
437  {
438  omegaup_s = 1.0;
439  omegaup_n = 0.0;
440  }
441  else
442  {
443  omegaup_s = 0.0;
444  omegaup_n = 1.0;
445  }
446 
447  // integration factor
448  RF factor = it->weight() * ig.geometry().integrationElement(it->position());
449 
450  // convection term
451  RF term1 = (omegaup_s*u_s + omegaup_n*u_n) * normalflux *factor;
452  for (size_type i=0; i<lfsv_s.size(); i++)
453  r_s.accumulate(lfsu_s,i,term1 * psi_s[i]);
454  for (size_type i=0; i<lfsv_n.size(); i++)
455  r_n.accumulate(lfsu_n,i,-term1 * psi_n[i]);
456 
457  // diffusion term
458  RF term2 = -(omega_s*(An_F_s*gradu_s) + omega_n*(An_F_n*gradu_n)) * factor;
459  for (size_type i=0; i<lfsv_s.size(); i++)
460  r_s.accumulate(lfsv_s,i,term2 * psi_s[i]);
461  for (size_type i=0; i<lfsv_n.size(); i++)
462  r_n.accumulate(lfsv_n,i,-term2 * psi_n[i]);
463 
464  // (non-)symmetric IP term
465  RF term3 = (u_s-u_n) * factor;
466  for (size_type i=0; i<lfsv_s.size(); i++)
467  r_s.accumulate(lfsv_s,i,term3 * theta * omega_s * (An_F_s*tgradpsi_s[i]));
468  for (size_type i=0; i<lfsv_n.size(); i++)
469  r_n.accumulate(lfsv_n,i,term3 * theta * omega_n * (An_F_n*tgradpsi_n[i]));
470 
471  // standard IP term integral
472  RF term4 = penalty_factor * (u_s-u_n) * factor;
473  for (size_type i=0; i<lfsv_s.size(); i++)
474  r_s.accumulate(lfsv_s,i,term4 * psi_s[i]);
475  for (size_type i=0; i<lfsv_n.size(); i++)
476  r_n.accumulate(lfsv_n,i,-term4 * psi_n[i]);
477  }
478  }
479 
480  template<typename IG, typename LFSU, typename X, typename LFSV, typename M>
481  void jacobian_skeleton (const IG& ig,
482  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
483  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
484  M& mat_ss, M& mat_sn,
485  M& mat_ns, M& mat_nn) const
486  {
487  // domain and range field type
488  typedef typename LFSV::Traits::FiniteElementType::
489  Traits::LocalBasisType::Traits::DomainFieldType DF;
490  typedef typename LFSV::Traits::FiniteElementType::
491  Traits::LocalBasisType::Traits::RangeFieldType RF;
492  typedef typename LFSV::Traits::FiniteElementType::
493  Traits::LocalBasisType::Traits::RangeType RangeType;
494  typedef typename LFSU::Traits::FiniteElementType::
495  Traits::LocalBasisType::Traits::JacobianType JacobianType;
496  typedef typename LFSV::Traits::SizeType size_type;
497 
498  // dimensions
499  const int dim = IG::dimension;
500  const int order = std::max(
501  std::max(lfsu_s.finiteElement().localBasis().order(),
502  lfsu_n.finiteElement().localBasis().order()),
503  std::max(lfsv_s.finiteElement().localBasis().order(),
504  lfsv_n.finiteElement().localBasis().order())
505  );
506  const int intorder = intorderadd+quadrature_factor*order;
507 
508  // evaluate permeability tensors
509  const Dune::FieldVector<DF,dim>&
510  inside_local = Dune::ReferenceElements<DF,dim>::general(ig.inside()->type()).position(0,0);
511  const Dune::FieldVector<DF,dim>&
512  outside_local = Dune::ReferenceElements<DF,dim>::general(ig.outside()->type()).position(0,0);
513  typename T::Traits::PermTensorType A_s, A_n;
514  A_s = param.A(*(ig.inside()),inside_local);
515  A_n = param.A(*(ig.outside()),outside_local);
516 
517  // face diameter; this should be revised for anisotropic meshes?
518  DF h_s, h_n;
519  DF hmax_s = 0., hmax_n = 0.;
520  element_size(ig.inside()->geometry(),h_s,hmax_s);
521  element_size(ig.outside()->geometry(),h_n,hmax_n);
522  RF h_F = std::min(h_s,h_n);
523  h_F = std::min(ig.inside()->geometry().volume(),ig.outside()->geometry().volume())/ig.geometry().volume(); // Houston!
524 
525  // select quadrature rule
526  Dune::GeometryType gtface = ig.geometryInInside().type();
527  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
528 
529  // transformation
530  typename IG::Entity::Geometry::JacobianInverseTransposed jac;
531 
532  // tensor times normal
533  const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
534  Dune::FieldVector<RF,dim> An_F_s;
535  A_s.mv(n_F,An_F_s);
536  Dune::FieldVector<RF,dim> An_F_n;
537  A_n.mv(n_F,An_F_n);
538 
539  // compute weights
540  RF omega_s;
541  RF omega_n;
542  RF harmonic_average(0.0);
544  {
545  RF delta_s = (An_F_s*n_F);
546  RF delta_n = (An_F_n*n_F);
547  omega_s = delta_n/(delta_s+delta_n+1e-20);
548  omega_n = delta_s/(delta_s+delta_n+1e-20);
549  harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1e-20);
550  }
551  else
552  {
553  omega_s = omega_n = 0.5;
554  harmonic_average = 1.0;
555  }
556 
557  // get polynomial degree
558  const int order_s = lfsu_s.finiteElement().localBasis().order();
559  const int order_n = lfsu_n.finiteElement().localBasis().order();
560  int degree = std::max( order_s, order_n );
561 
562  // penalty factor
563  RF penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
564 
565  // loop over quadrature points
566  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
567  {
568  // exact normal
569  const Dune::FieldVector<DF,dim> n_F_local = ig.unitOuterNormal(it->position());
570 
571  // position of quadrature point in local coordinates of elements
572  Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(it->position());
573  Dune::FieldVector<DF,dim> iplocal_n = ig.geometryInOutside().global(it->position());
574 
575  // evaluate basis functions
576 #if USECACHE==0
577  std::vector<RangeType> phi_s(lfsu_s.size());
578  lfsu_s.finiteElement().localBasis().evaluateFunction(iplocal_s,phi_s);
579  std::vector<RangeType> phi_n(lfsu_n.size());
580  lfsu_n.finiteElement().localBasis().evaluateFunction(iplocal_n,phi_n);
581 #else
582  const std::vector<RangeType>& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
583  const std::vector<RangeType>& phi_n = cache[order_n].evaluateFunction(iplocal_n,lfsu_n.finiteElement().localBasis());
584 #endif
585 
586  // evaluate gradient of basis functions (we assume Galerkin method lfsu=lfsv)
587 #if USECACHE==0
588  std::vector<JacobianType> gradphi_s(lfsu_s.size());
589  lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
590  std::vector<JacobianType> gradphi_n(lfsu_n.size());
591  lfsu_n.finiteElement().localBasis().evaluateJacobian(iplocal_n,gradphi_n);
592 #else
593  const std::vector<JacobianType>& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
594  const std::vector<JacobianType>& gradphi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsu_n.finiteElement().localBasis());
595 #endif
596 
597  // transform gradients of shape functions to real element
598  jac = ig.inside()->geometry().jacobianInverseTransposed(iplocal_s);
599  std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
600  for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
601  jac = ig.outside()->geometry().jacobianInverseTransposed(iplocal_n);
602  std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
603  for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
604 
605  // evaluate velocity field and upwinding, assume H(div) velocity field => may choose any side
606  typename T::Traits::RangeType b = param.b(*(ig.inside()),iplocal_s);
607  RF normalflux = b*n_F_local;
608  RF omegaup_s, omegaup_n;
609  if (normalflux>=0.0)
610  {
611  omegaup_s = 1.0;
612  omegaup_n = 0.0;
613  }
614  else
615  {
616  omegaup_s = 0.0;
617  omegaup_n = 1.0;
618  }
619 
620  // integration factor
621  RF factor = it->weight() * ig.geometry().integrationElement(it->position());
622  RF ipfactor = penalty_factor * factor;
623 
624  // do all terms in the order: I convection, II diffusion, III consistency, IV ip
625  for (size_type j=0; j<lfsu_s.size(); j++) {
626  RF temp1 = -(An_F_s*tgradphi_s[j])*omega_s*factor;
627  for (size_type i=0; i<lfsu_s.size(); i++) {
628  mat_ss.accumulate(lfsu_s,i,lfsu_s,j,omegaup_s * phi_s[j] * normalflux *factor * phi_s[i]);
629  mat_ss.accumulate(lfsu_s,i,lfsu_s,j,temp1 * phi_s[i]);
630  mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * factor * theta * omega_s * (An_F_s*tgradphi_s[i]));
631  mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * ipfactor * phi_s[i]);
632  }
633  }
634  for (size_type j=0; j<lfsu_n.size(); j++) {
635  RF temp1 = -(An_F_n*tgradphi_n[j])*omega_n*factor;
636  for (size_type i=0; i<lfsu_s.size(); i++) {
637  mat_sn.accumulate(lfsu_s,i,lfsu_n,j,omegaup_n * phi_n[j] * normalflux *factor * phi_s[i]);
638  mat_sn.accumulate(lfsu_s,i,lfsu_n,j,temp1 * phi_s[i]);
639  mat_sn.accumulate(lfsu_s,i,lfsu_n,j,-phi_n[j] * factor * theta * omega_s * (An_F_s*tgradphi_s[i]));
640  mat_sn.accumulate(lfsu_s,i,lfsu_n,j,-phi_n[j] * ipfactor * phi_s[i]);
641  }
642  }
643  for (size_type j=0; j<lfsu_s.size(); j++) {
644  RF temp1 = -(An_F_s*tgradphi_s[j])*omega_s*factor;
645  for (size_type i=0; i<lfsu_n.size(); i++) {
646  mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-omegaup_s * phi_s[j] * normalflux *factor * phi_n[i]);
647  mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-temp1 * phi_n[i]);
648  mat_ns.accumulate(lfsu_n,i,lfsu_s,j,phi_s[j] * factor * theta * omega_n * (An_F_n*tgradphi_n[i]));
649  mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-phi_s[j] * ipfactor * phi_n[i]);
650  }
651  }
652  for (size_type j=0; j<lfsu_n.size(); j++) {
653  RF temp1 = -(An_F_n*tgradphi_n[j])*omega_n*factor;
654  for (size_type i=0; i<lfsu_n.size(); i++) {
655  mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-omegaup_n * phi_n[j] * normalflux *factor * phi_n[i]);
656  mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-temp1 * phi_n[i]);
657  mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-phi_n[j] * factor * theta * omega_n * (An_F_n*tgradphi_n[i]));
658  mat_nn.accumulate(lfsu_n,i,lfsu_n,j,phi_n[j] * ipfactor * phi_n[i]);
659  }
660  }
661  }
662  }
663 
664  // boundary integral depending on test and ansatz functions
665  // We put the Dirchlet evaluation also in the alpha term to save some geometry evaluations
666  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
667  void alpha_boundary (const IG& ig,
668  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
669  R& r_s) const
670  {
671  // domain and range field type
672  typedef typename LFSV::Traits::FiniteElementType::
673  Traits::LocalBasisType::Traits::DomainFieldType DF;
674  typedef typename LFSV::Traits::FiniteElementType::
675  Traits::LocalBasisType::Traits::RangeFieldType RF;
676  typedef typename LFSV::Traits::FiniteElementType::
677  Traits::LocalBasisType::Traits::RangeType RangeType;
678  typedef typename LFSU::Traits::FiniteElementType::
679  Traits::LocalBasisType::Traits::JacobianType JacobianType;
680  typedef typename LFSV::Traits::SizeType size_type;
681 
682  // dimensions
683  const int dim = IG::dimension;
684  const int order = std::max(
685  lfsu_s.finiteElement().localBasis().order(),
686  lfsv_s.finiteElement().localBasis().order()
687  );
688  const int intorder = intorderadd+quadrature_factor*order;
689 
690  // evaluate permeability tensors
691  const Dune::FieldVector<DF,dim>&
692  inside_local = Dune::ReferenceElements<DF,dim>::general(ig.inside()->type()).position(0,0);
693  typename T::Traits::PermTensorType A_s;
694  A_s = param.A(*(ig.inside()),inside_local);
695 
696  // select quadrature rule
697  Dune::GeometryType gtface = ig.geometryInInside().type();
698  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
699 
700  // face diameter
701  DF h_s;
702  DF hmax_s = 0.;
703  element_size(ig.inside()->geometry(),h_s,hmax_s);
704  RF h_F = h_s;
705  h_F = ig.inside()->geometry().volume()/ig.geometry().volume(); // Houston!
706 
707  // transformation
708  Dune::FieldMatrix<DF,dim,dim> jac;
709 
710  // compute weights
711  const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
712  Dune::FieldVector<RF,dim> An_F_s;
713  A_s.mv(n_F,An_F_s);
714  RF harmonic_average;
716  harmonic_average = An_F_s*n_F;
717  else
718  harmonic_average = 1.0;
719 
720  // get polynomial degree
721  const int order_s = lfsu_s.finiteElement().localBasis().order();
722  int degree = order_s;
723 
724  // penalty factor
725  RF penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
726 
727  // loop over quadrature points
728  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
729  {
730  BCType bctype = param.bctype(ig.intersection(),it->position());
731 
732  // position of quadrature point in local coordinates of elements
733  Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(it->position());
734 
735  // local normal
736  const Dune::FieldVector<DF,dim> n_F_local = ig.unitOuterNormal(it->position());
737 
738  // evaluate basis functions
739 #if USECACHE==0
740  std::vector<RangeType> phi_s(lfsu_s.size());
741  lfsu_s.finiteElement().localBasis().evaluateFunction(iplocal_s,phi_s);
742  std::vector<RangeType> psi_s(lfsv_s.size());
743  lfsv_s.finiteElement().localBasis().evaluateFunction(iplocal_s,psi_s);
744 #else
745  const std::vector<RangeType>& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
746  const std::vector<RangeType>& psi_s = cache[order_s].evaluateFunction(iplocal_s,lfsv_s.finiteElement().localBasis());
747 #endif
748 
749  // integration factor
750  RF factor = it->weight() * ig.geometry().integrationElement(it->position());
751 
753  {
754  // evaluate flux boundary condition
755  RF j = param.j(ig.intersection(),it->position());
756 
757  // integrate
758  for (size_type i=0; i<lfsv_s.size(); i++)
759  r_s.accumulate(lfsv_s,i,j * psi_s[i] * factor);
760 
761  continue;
762  }
763 
764  // evaluate u
765  RF u_s=0.0;
766  for (size_type i=0; i<lfsu_s.size(); i++)
767  u_s += x_s(lfsu_s,i)*phi_s[i];
768 
769  // evaluate velocity field and upwinding, assume H(div) velocity field => choose any side
770  typename T::Traits::RangeType b = param.b(*(ig.inside()),iplocal_s);
771  RF normalflux = b*n_F_local;
772 
774  {
775  if (normalflux<-1e-30)
776  DUNE_THROW(Dune::Exception,
777  "Outflow boundary condition on inflow! [b("
778  << ig.geometry().global(it->position()) << ") = "
779  << b << ")");
780 
781  // convection term
782  RF term1 = u_s * normalflux *factor;
783  for (size_type i=0; i<lfsv_s.size(); i++)
784  r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
785 
786  // evaluate flux boundary condition
787  RF o = param.o(ig.intersection(),it->position());
788 
789  // integrate
790  for (size_type i=0; i<lfsv_s.size(); i++)
791  r_s.accumulate(lfsv_s,i,o * psi_s[i] * factor);
792 
793  continue;
794  }
795 
796  // evaluate gradient of basis functions (we assume Galerkin method lfsu=lfsv)
798 #if USECACHE==0
799  std::vector<JacobianType> gradphi_s(lfsu_s.size());
800  lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
801  std::vector<JacobianType> gradpsi_s(lfsv_s.size());
802  lfsv_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradpsi_s);
803 #else
804  const std::vector<JacobianType>& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
805  const std::vector<JacobianType>& gradpsi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsv_s.finiteElement().localBasis());
806 #endif
807 
808  // transform gradients of shape functions to real element
809  jac = ig.inside()->geometry().jacobianInverseTransposed(iplocal_s);
810  std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
811  for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
812  std::vector<Dune::FieldVector<RF,dim> > tgradpsi_s(lfsv_s.size());
813  for (size_type i=0; i<lfsv_s.size(); i++) jac.mv(gradpsi_s[i][0],tgradpsi_s[i]);
814 
815  // compute gradient of u
816  Dune::FieldVector<RF,dim> gradu_s(0.0);
817  for (size_type i=0; i<lfsu_s.size(); i++)
818  gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
819 
820  // evaluate Dirichlet boundary condition
821  RF g = param.g(*(ig.inside()),iplocal_s);
822 
823  // upwind
824  RF omegaup_s, omegaup_n;
825  if (normalflux>=0.0)
826  {
827  omegaup_s = 1.0;
828  omegaup_n = 0.0;
829  }
830  else
831  {
832  omegaup_s = 0.0;
833  omegaup_n = 1.0;
834  }
835 
836  // convection term
837  RF term1 = (omegaup_s*u_s + omegaup_n*g) * normalflux *factor;
838  for (size_type i=0; i<lfsv_s.size(); i++)
839  r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
840 
841  // diffusion term
842  RF term2 = (An_F_s*gradu_s) * factor;
843  for (size_type i=0; i<lfsv_s.size(); i++)
844  r_s.accumulate(lfsv_s,i,-term2 * psi_s[i]);
845 
846  // (non-)symmetric IP term
847  RF term3 = (u_s-g) * factor;
848  for (size_type i=0; i<lfsv_s.size(); i++)
849  r_s.accumulate(lfsv_s,i,term3 * theta * (An_F_s*tgradpsi_s[i]));
850 
851  // standard IP term
852  RF term4 = penalty_factor * (u_s-g) * factor;
853  for (size_type i=0; i<lfsv_s.size(); i++)
854  r_s.accumulate(lfsv_s,i,term4 * psi_s[i]);
855  }
856  }
857 
858  template<typename IG, typename LFSU, typename X, typename LFSV, typename M>
859  void jacobian_boundary (const IG& ig,
860  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
861  M& mat_ss) const
862  {
863  // domain and range field type
864  typedef typename LFSV::Traits::FiniteElementType::
865  Traits::LocalBasisType::Traits::DomainFieldType DF;
866  typedef typename LFSV::Traits::FiniteElementType::
867  Traits::LocalBasisType::Traits::RangeFieldType RF;
868  typedef typename LFSV::Traits::FiniteElementType::
869  Traits::LocalBasisType::Traits::RangeType RangeType;
870  typedef typename LFSU::Traits::FiniteElementType::
871  Traits::LocalBasisType::Traits::JacobianType JacobianType;
872  typedef typename LFSV::Traits::SizeType size_type;
873 
874  // dimensions
875  const int dim = IG::dimension;
876  const int order = std::max(
877  lfsu_s.finiteElement().localBasis().order(),
878  lfsv_s.finiteElement().localBasis().order()
879  );
880  const int intorder = intorderadd+quadrature_factor*order;
881 
882  // evaluate permeability tensors
883  const Dune::FieldVector<DF,dim>&
884  inside_local = Dune::ReferenceElements<DF,dim>::general(ig.inside()->type()).position(0,0);
885  typename T::Traits::PermTensorType A_s;
886  A_s = param.A(*(ig.inside()),inside_local);
887 
888  // select quadrature rule
889  Dune::GeometryType gtface = ig.geometryInInside().type();
890  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
891 
892  // face diameter
893  DF h_s;
894  DF hmax_s = 0.;
895  element_size(ig.inside()->geometry(),h_s,hmax_s);
896  RF h_F = h_s;
897  h_F = ig.inside()->geometry().volume()/ig.geometry().volume(); // Houston!
898 
899  // transformation
900  Dune::FieldMatrix<DF,dim,dim> jac;
901 
902  // compute weights
903  const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
904  Dune::FieldVector<RF,dim> An_F_s;
905  A_s.mv(n_F,An_F_s);
906  RF harmonic_average;
908  harmonic_average = An_F_s*n_F;
909  else
910  harmonic_average = 1.0;
911 
912  // get polynomial degree
913  const int order_s = lfsu_s.finiteElement().localBasis().order();
914  int degree = order_s;
915 
916  // penalty factor
917  RF penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
918 
919  // Neumann boundary makes no contribution to boundary
920  //if (bctype == ConvectionDiffusionBoundaryConditions::Neumann) return;
921 
922  // loop over quadrature points
923  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
924  {
925  BCType bctype = param.bctype(ig.intersection(),it->position());
926  if (bctype == ConvectionDiffusionBoundaryConditions::Neumann) continue;
927 
928  // position of quadrature point in local coordinates of elements
929  Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(it->position());
930 
931  // local normal
932  const Dune::FieldVector<DF,dim> n_F_local = ig.unitOuterNormal(it->position());
933 
934  // evaluate basis functions
935 #if USECACHE==0
936  std::vector<RangeType> phi_s(lfsu_s.size());
937  lfsu_s.finiteElement().localBasis().evaluateFunction(iplocal_s,phi_s);
938 #else
939  const std::vector<RangeType>& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
940 #endif
941 
942  // integration factor
943  RF factor = it->weight() * ig.geometry().integrationElement(it->position());
944 
945  // evaluate velocity field and upwinding, assume H(div) velocity field => choose any side
946  typename T::Traits::RangeType b = param.b(*(ig.inside()),iplocal_s);
947  RF normalflux = b*n_F_local;
948 
950  {
951  if (normalflux<-1e-30)
952  DUNE_THROW(Dune::Exception,
953  "Outflow boundary condition on inflow! [b("
954  << ig.geometry().global(it->position()) << ") = "
955  << b << ")" << n_F_local << " " << normalflux);
956 
957  // convection term
958  for (size_type j=0; j<lfsu_s.size(); j++)
959  for (size_type i=0; i<lfsu_s.size(); i++)
960  mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * normalflux * factor * phi_s[i]);
961 
962  continue;
963  }
964 
965  // evaluate gradient of basis functions (we assume Galerkin method lfsu=lfsv)
966 #if USECACHE==0
967  std::vector<JacobianType> gradphi_s(lfsu_s.size());
968  lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
969 #else
970  const std::vector<JacobianType>& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
971 #endif
972 
973  // transform gradients of shape functions to real element
974  jac = ig.inside()->geometry().jacobianInverseTransposed(iplocal_s);
975  std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
976  for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
977 
978  // upwind
979  RF omegaup_s, omegaup_n;
980  if (normalflux>=0.0)
981  {
982  omegaup_s = 1.0;
983  omegaup_n = 0.0;
984  }
985  else
986  {
987  omegaup_s = 0.0;
988  omegaup_n = 1.0;
989  }
990 
991  // convection term
992  for (size_type j=0; j<lfsu_s.size(); j++)
993  for (size_type i=0; i<lfsu_s.size(); i++)
994  mat_ss.accumulate(lfsu_s,i,lfsu_s,j,omegaup_s * phi_s[j] * normalflux * factor * phi_s[i]);
995 
996  // diffusion term
997  for (size_type j=0; j<lfsu_s.size(); j++)
998  for (size_type i=0; i<lfsu_s.size(); i++)
999  mat_ss.accumulate(lfsu_s,i,lfsu_s,j,-(An_F_s*tgradphi_s[j]) * factor * phi_s[i]);
1000 
1001  // (non-)symmetric IP term
1002  for (size_type j=0; j<lfsu_s.size(); j++)
1003  for (size_type i=0; i<lfsu_s.size(); i++)
1004  mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * factor * theta * (An_F_s*tgradphi_s[i]));
1005 
1006  // standard IP term
1007  for (size_type j=0; j<lfsu_s.size(); j++)
1008  for (size_type i=0; i<lfsu_s.size(); i++)
1009  mat_ss.accumulate(lfsu_s,i,lfsu_s,j,penalty_factor * phi_s[j] * phi_s[i] * factor);
1010  }
1011  }
1012 
1013  // volume integral depending only on test functions
1014  template<typename EG, typename LFSV, typename R>
1015  void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
1016  {
1017  // domain and range field type
1018  typedef typename LFSV::Traits::FiniteElementType::
1019  Traits::LocalBasisType::Traits::DomainFieldType DF;
1020  typedef typename LFSV::Traits::FiniteElementType::
1021  Traits::LocalBasisType::Traits::RangeFieldType RF;
1022  typedef typename LFSV::Traits::FiniteElementType::
1023  Traits::LocalBasisType::Traits::RangeType RangeType;
1024  typedef typename LFSV::Traits::SizeType size_type;
1025 
1026  // dimensions
1027  const int dim = EG::Geometry::dimension;
1028  const int order = lfsv.finiteElement().localBasis().order();
1029  const int intorder = intorderadd + 2 * order;
1030 
1031  // select quadrature rule
1032  Dune::GeometryType gt = eg.geometry().type();
1033  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
1034 
1035  // loop over quadrature points
1036  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
1037  {
1038  // evaluate shape functions
1039 #if USECACHE==0
1040  std::vector<RangeType> phi(lfsv.size());
1041  lfsv.finiteElement().localBasis().evaluateFunction(it->position(),phi);
1042 #else
1043  const std::vector<RangeType>& phi = cache[order].evaluateFunction(it->position(),lfsv.finiteElement().localBasis());
1044 #endif
1045 
1046  // evaluate right hand side parameter function
1047  Real f;
1048  f = param.f(eg.entity(),it->position());
1049 
1050  // integrate f
1051  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
1052  for (size_type i=0; i<lfsv.size(); i++)
1053  r.accumulate(lfsv,i,-f*phi[i]*factor);
1054  }
1055  }
1056 
1058  void setTime (double t)
1059  {
1060  param.setTime(t);
1061  }
1062 
1063  private:
1064  T& param; // two phase parameter class
1067  Real alpha, beta;
1068  int intorderadd;
1069  int quadrature_factor;
1070  Real theta;
1071  typedef typename FiniteElementMap::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
1072 
1074 
1075  // In theory it is possible that one and the same local operator is
1076  // called first with a finite element of one type and later with a
1077  // finite element of another type. Since finite elements of different
1078  // type will usually produce different results for the same local
1079  // coordinate they cannot share a cache. Here we use a vector of caches
1080  // to allow for different orders of the shape functions, which should be
1081  // enough to support p-adaptivity. (Another likely candidate would be
1082  // differing geometry types, i.e. hybrid meshes.)
1083 
1084  std::vector<Cache> cache;
1085 
1086  template<class GEO>
1087  void element_size (const GEO& geo, typename GEO::ctype& hmin, typename GEO::ctype hmax) const
1088  {
1089  typedef typename GEO::ctype DF;
1090  hmin = 1.0E100;
1091  hmax = -1.0E00;
1092  const int dim = GEO::coorddimension;
1093  if (dim==1)
1094  {
1095  Dune::FieldVector<DF,dim> x = geo.corner(0);
1096  x -= geo.corner(1);
1097  hmin = hmax = x.two_norm();
1098  return;
1099  }
1100  else
1101  {
1102  Dune::GeometryType gt = geo.type();
1103  for (int i=0; i<Dune::ReferenceElements<DF,dim>::general(gt).size(dim-1); i++)
1104  {
1105  Dune::FieldVector<DF,dim> x = geo.corner(Dune::ReferenceElements<DF,dim>::general(gt).subEntity(i,dim-1,0,dim));
1106  x -= geo.corner(Dune::ReferenceElements<DF,dim>::general(gt).subEntity(i,dim-1,1,dim));
1107  hmin = std::min(hmin,x.two_norm());
1108  hmax = std::max(hmax,x.two_norm());
1109  }
1110  return;
1111  }
1112  }
1113  };
1114  }
1115 }
1116 #endif
Definition: convectiondiffusionparameter.hh:68
Definition: convectiondiffusiondg.hh:35
sparsity pattern generator
Definition: pattern.hh:30
Implement jacobian_apply_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:445
Definition: convectiondiffusiondg.hh:74
Default flags for all local operators.
Definition: flags.hh:18
Definition: convectiondiffusiondg.hh:80
Type
Definition: convectiondiffusiondg.hh:40
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538
void setTime(double t)
set time in parameter class
Definition: convectiondiffusiondg.hh:1058
void alpha_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, R &r_s, R &r_n) const
Definition: convectiondiffusiondg.hh:276
Definition: convectiondiffusiondg.hh:40
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: convectiondiffusiondg.hh:1015
Definition: convectiondiffusiondg.hh:78
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusiondg.hh:103
Definition: convectiondiffusiondg.hh:58
const IG & ig
Definition: common/constraints.hh:146
Definition: convectiondiffusiondg.hh:40
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusiondg.hh:198
void jacobian_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, M &mat_ss) const
Definition: convectiondiffusiondg.hh:859
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: convectiondiffusiondg.hh:667
Definition: convectiondiffusiondg.hh:38
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:16
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
const F & f
Definition: common/constraints.hh:145
Definition: convectiondiffusiondg.hh:75
void jacobian_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, M &mat_ss, M &mat_sn, M &mat_ns, M &mat_nn) const
Definition: convectiondiffusiondg.hh:481
sparsity pattern generator
Definition: pattern.hh:14
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
ConvectionDiffusionDG(T &param_, ConvectionDiffusionDGMethod::Type method_=ConvectionDiffusionDGMethod::NIPG, ConvectionDiffusionDGWeights::Type weights_=ConvectionDiffusionDGWeights::weightsOff, Real alpha_=0.0, int intorderadd_=0)
constructor: pass parameter object
Definition: convectiondiffusiondg.hh:84
const EG & eg
Definition: common/constraints.hh:277
Definition: convectiondiffusiondg.hh:33
Definition: convectiondiffusiondg.hh:79
Type
Definition: convectiondiffusiondg.hh:35
Definition: convectiondiffusiondg.hh:81
Definition: convectiondiffusionparameter.hh:68
const E & e
Definition: interpolate.hh:172
Type
Definition: convectiondiffusionparameter.hh:68
Definition: convectiondiffusionparameter.hh:68
Definition: convectiondiffusiondg.hh:35