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