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