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