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