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