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