dune-pdelab  2.4.1
diffusiondg.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_DIFFUSIONDG_HH
3 #define DUNE_PDELAB_DIFFUSIONDG_HH
4 #warning This file is deprecated and will be removed after the Dune-PDELab 2.4 release! Use the ConvectionDiffusionDG local operator from dune/pdelab/localoperator/convectiondiffusiondg.hh instead!
5 
6 #include <dune/common/exceptions.hh>
7 #include <dune/common/fvector.hh>
8 #include <dune/geometry/quadraturerules.hh>
9 #include <dune/geometry/referenceelements.hh>
10 
12 
13 #include "defaultimp.hh"
14 #include "pattern.hh"
15 #include "flags.hh"
16 #include "diffusionparam.hh"
17 
18 namespace Dune {
19  namespace PDELab {
20 
21  // a local operator for solving the diffusion equation
22  // - div (K grad u) = f in \Omega,
23  // u = g on \Gamma_D
24  // (- K grad u) * nu = j on \Gamma_N
25  // discontinuous Galerkin method (SIPG, NIPG, OBB)
26  //
27  // @tparam K grid function for permeability tensor
28  // @tparam F grid function for rhs
29  // @tparam B boundary type function
30  // @tparam G grid function for Dirichlet boundary conditions
31  // @tparam J grid function for Neumann boundary conditions
32  template<typename K, typename F, typename B, typename G, typename J>
33  class DiffusionDG :
36 // #define JacobianBasedAlphaX
37 // #define NumericalJacobianX
38 #ifdef JacobianBasedAlphaX
39  ,public JacobianBasedAlphaVolume<DiffusionDG<K, F, B, G, J> >
40  ,public JacobianBasedAlphaSkeleton<DiffusionDG<K, F, B, G, J> >
41  ,public JacobianBasedAlphaBoundary<DiffusionDG<K, F, B, G, J> >
42 #endif
43 #ifdef NumericalJacobianX
44  #ifdef JacobianBasedAlphaX
45  #error You have provide either the alpha_* or the jacobian_* methods...
46  #endif
47  ,public NumericalJacobianVolume<DiffusionDG<K, F, B, G, J> >
48  ,public NumericalJacobianSkeleton<DiffusionDG<K, F, B, G, J> >
49  ,public NumericalJacobianBoundary<DiffusionDG<K, F, B, G, J> >
50 #endif
51  {
52  public:
53  // pattern assembly flags
54  enum { doPatternVolume = true };
55  enum { doPatternSkeleton = true };
56 
57  // residual assembly flags
58  enum { doAlphaVolume = true };
59  enum { doAlphaSkeleton = true };
60  enum { doAlphaBoundary = true };
61  enum { doLambdaVolume = true };
62  enum { doLambdaSkeleton = false };
63  enum { doLambdaBoundary = true };
64 
65  DUNE_DEPRECATED_MSG("Deprecated in Dune-PDELab 2.4, use the local operator ConvectionDiffusionDG instead!")
66  DiffusionDG (const K& k_, const F& f_, const B& bctype_, const G& g_, const J& j_, int dg_method, int _superintegration_order = 0) :
67  k(k_), f(f_), bctype(bctype_), g(g_), j(j_), superintegration_order(_superintegration_order)
68  {
69 
70  // OBB
71  if (dg_method == 0)
72  {
73  epsilon = 1.0;
74  sigma = 0.0;
75  beta = 0.0;
76  }
77  // NIPG
78  else if (dg_method == 1)
79  {
80  epsilon = 1.0;
81  sigma = 4.5; // should be < 5 for stability reasons
82  beta = 2.0 - 0.5*F::Traits::dimDomain; // 2D => 1, 3D => 0.5
83  }
84  // SIPG
85  else if (dg_method == 2)
86  {
87  epsilon = -1.0;
88  sigma = 4.5; // should be < 5 for stability reasons
89  beta = 2.0 - 0.5*F::Traits::dimDomain; // 2D => 1, 3D => 0.5
90  }
91  }
92 
93 #ifndef JacobianBasedAlphaX
94  // volume integral depending on test and ansatz functions
95  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
96  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
97  {
98  // domain and range field type
99  typedef typename LFSU::Traits::FiniteElementType::
100  Traits::LocalBasisType::Traits::DomainFieldType DF;
101  typedef typename LFSU::Traits::FiniteElementType::
102  Traits::LocalBasisType::Traits::RangeFieldType RF;
103  typedef typename LFSU::Traits::FiniteElementType::
104  Traits::LocalBasisType::Traits::JacobianType JacobianType;
105 
106  // dimensionslocal
107  const int dim = EG::Geometry::dimension;
108 
109  // select quadrature rule
110  Dune::GeometryType gt = eg.geometry().type();
111  const int qorder = std::max ( 2 * ( (int)lfsu.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
112  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
113 
114  // evaluate diffusion tensor at cell center, assume it is constant over elements
115  typename K::Traits::RangeType tensor(0.0);
116  Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
117  k.evaluate(eg.entity(),localcenter,tensor);
118 
119  // loop over quadrature points
120  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
121  {
122  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
123  std::vector<JacobianType> js(lfsu.size());
124  lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
125 
126  // transform gradient to real element
127  const typename EG::Geometry::JacobianInverseTransposed jac =
128  eg.geometry().jacobianInverseTransposed(it->position());
129  std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
130  for (size_t i=0; i<lfsu.size(); i++)
131  {
132  gradphi[i] = 0.0;
133  jac.umv(js[i][0],gradphi[i]);
134  }
135 
136  // compute gradient of u
137  Dune::FieldVector<RF,dim> gradu(0.0);
138  for (size_t i=0; i<lfsu.size(); i++)
139  gradu.axpy(x(lfsu,i),gradphi[i]);
140 
141  // compute K * gradient of u
142  Dune::FieldVector<RF,dim> Kgradu(0.0);
143  tensor.umv(gradu,Kgradu);
144 
145  // integrate (K grad u)*grad phi_i
146  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
147  for (size_t i=0; i<lfsu.size(); i++)
148  {
149  r.accumulate( lfsv, i, Kgradu*gradphi[i] * factor ) ;
150  }
151  }
152  }
153 
154  // skeleton integral depending on test and ansatz functions
155  // each face is only visited ONCE!
156  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
157  void alpha_skeleton (const IG& ig,
158  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
159  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
160  R& r_s, R& r_n) const
161  {
162  // domain and range field type
163  typedef typename LFSU::Traits::FiniteElementType::
164  Traits::LocalBasisType::Traits::DomainFieldType DF;
165  typedef typename LFSU::Traits::FiniteElementType::
166  Traits::LocalBasisType::Traits::RangeFieldType RF;
167  typedef typename LFSV::Traits::FiniteElementType::
168  Traits::LocalBasisType::Traits::RangeType RangeType;
169  typedef typename LFSV::Traits::FiniteElementType::
170  Traits::LocalBasisType::Traits::JacobianType JacobianType;
171 
172  // dimensions
173  const int dim = IG::dimension;
174  const int dimw = IG::dimensionworld;
175 
176  // select quadrature rule
177  Dune::GeometryType gtface = ig.geometryInInside().type();
178  const int qorder = std::max( 0, std::max(
179  2 * ( (int)lfsu_s.finiteElement().localBasis().order() - 1 ),
180  2 * ( (int)lfsu_n.finiteElement().localBasis().order() - 1 ))) + superintegration_order;
181  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
182 
183  // normal of center in face's reference element
184  const Dune::FieldVector<DF,IG::dimension-1>& face_center =
185  Dune::ReferenceElements<DF,IG::dimension-1>::
186  general(ig.geometry().type()).position(0,0);
187  const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
188 
189  // evaluate diffusion tensor at elements' centers, assume they are constant over elements
190  const Dune::FieldVector<DF,IG::dimension>&
191  inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
192  const Dune::FieldVector<DF,IG::dimension>&
193  outside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.outside()->type()).position(0,0);
194  typename K::Traits::RangeType permeability_s(0.0);
195  typename K::Traits::RangeType permeability_n(0.0);
196  k.evaluate(*(ig.inside()),inside_local,permeability_s);
197  k.evaluate(*(ig.outside()),outside_local,permeability_n);
198 
199  // penalty weight for NIPG / SIPG
200  RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
201 
202  // loop over quadrature points
203  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
204  {
205  // position of quadrature point in local coordinates of element
206  Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(it->position());
207  Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(it->position());
208 
209  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
210  std::vector<JacobianType> js_s(lfsv_s.size());
211  lfsv_s.finiteElement().localBasis().evaluateJacobian(local_s,js_s);
212  std::vector<JacobianType> js_n(lfsv_n.size());
213  lfsv_n.finiteElement().localBasis().evaluateJacobian(local_n,js_n);
214 
215  // transform gradient to real element
216  typename IG::Entity::Geometry::JacobianInverseTransposed jac_s;
217  jac_s = ig.inside()->geometry().jacobianInverseTransposed(local_s);
218  std::vector<Dune::FieldVector<RF,dim> > gradphi_s(lfsv_s.size());
219  for (size_t i=0; i<lfsv_s.size(); i++)
220  {
221  gradphi_s[i] = 0.0;
222  jac_s.umv(js_s[i][0],gradphi_s[i]);
223  }
224  typename IG::Entity::Geometry::JacobianInverseTransposed jac_n;
225  jac_n = ig.outside()->geometry().jacobianInverseTransposed(local_n);
226  std::vector<Dune::FieldVector<RF,dim> > gradphi_n(lfsv_n.size());
227  for (size_t i=0; i<lfsv_n.size(); i++)
228  {
229  gradphi_n[i] = 0.0;
230  jac_n.umv(js_n[i][0],gradphi_n[i]);
231  }
232 
233  // evaluate test shape functions
234  std::vector<RangeType> phi_s(lfsv_s.size());
235  lfsv_s.finiteElement().localBasis().evaluateFunction(local_s,phi_s);
236  std::vector<RangeType> phi_n(lfsv_n.size());
237  lfsv_n.finiteElement().localBasis().evaluateFunction(local_n,phi_n);
238 
239  // compute gradient of u
240  Dune::FieldVector<RF,dim> gradu_s(0.0);
241  for (size_t i=0; i<lfsu_s.size(); i++)
242  {
243  gradu_s.axpy(x_s(lfsu_s,i),gradphi_s[i]);
244  }
245  Dune::FieldVector<RF,dim> gradu_n(0.0);
246  for (size_t i=0; i<lfsu_n.size(); i++)
247  {
248  gradu_n.axpy(x_n(lfsu_n,i),gradphi_n[i]);
249  }
250 
251  // compute K * gradient of u
252  Dune::FieldVector<RF,dim> kgradu_s(0.0);
253  permeability_s.umv(gradu_s,kgradu_s);
254  Dune::FieldVector<RF,dim> kgradu_n(0.0);
255  permeability_n.umv(gradu_n,kgradu_n);
256 
257  // evaluate u in both elements self and neighbor
258  RF u_s = 0.0;
259  for (size_t i=0; i<lfsu_s.size(); i++)
260  {
261  u_s += x_s(lfsu_s,i)*phi_s[i];
262  }
263  RF u_n = 0.0;
264  for (size_t i=0; i<lfsu_n.size(); i++)
265  {
266  u_n += x_n(lfsu_n,i)*phi_n[i];
267  }
268 
269  // jump and average for u
270  RF u_jump = u_s - u_n;
271 
272  // average on intersection of K * grad u * normal
273  RF kgradunormal_average = (kgradu_s + kgradu_n)*normal * 0.5;
274 
275  // average on intersection of K * grad v * normal
276  std::vector<Dune::FieldVector<RF,dim> > kgradphi_s(lfsu_s.size());
277  std::vector<Dune::FieldVector<RF,dim> > kgradphi_n(lfsu_n.size());
278  for (size_t i=0; i<lfsu_s.size(); i++)
279  {
280  permeability_s.mv(gradphi_s[i],kgradphi_s[i]);
281  }
282  for (size_t i=0; i<lfsu_n.size(); i++)
283  {
284  permeability_n.mv(gradphi_n[i],kgradphi_n[i]);
285  }
286 
287  // integrate what needed
288  RF factor = it->weight()*ig.geometry().integrationElement(it->position());
289  for (unsigned int i=0; i<lfsv_s.size(); i++)
290  {
291  // NIPG / SIPG penalty term: sigma/|gamma|^beta * [u]*[v]
292  r_s.accumulate( lfsv_s, i, penalty_weight * u_jump*phi_s[i]*factor );
293  // epsilon * <Kgradv*my>[u] - <Kgradu*my>[v]
294  r_s.accumulate( lfsv_s, i, epsilon*(kgradphi_s[i]*normal)*0.5*u_jump*factor );
295  r_s.accumulate( lfsv_s, i, - phi_s[i]*kgradunormal_average*factor );
296  }
297  for (unsigned int i=0; i<lfsv_n.size(); i++)
298  {
299  // NIPG / SIPG penalty term: sigma/|gamma|^beta * [u]*[v]
300  r_n.accumulate( lfsv_s, i, penalty_weight * u_jump*(-phi_n[i])*factor );
301  // epsilon * <Kgradv*my>[u] - [v]<Kgradu*my>
302  r_n.accumulate( lfsv_s, i, epsilon*(kgradphi_n[i]*normal)*0.5*u_jump*factor );
303  r_n.accumulate( lfsv_s, i, phi_n[i] * kgradunormal_average * factor );
304  }
305  }
306  }
307 
308  // boundary integral depending on test and ansatz functions
309  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
310  void alpha_boundary (const IG& ig, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
311  {
312  // domain and range field type
313  typedef typename LFSU::Traits::FiniteElementType::
314  Traits::LocalBasisType::Traits::DomainFieldType DF;
315  typedef typename LFSU::Traits::FiniteElementType::
316  Traits::LocalBasisType::Traits::RangeFieldType RF;
317  typedef typename LFSV::Traits::FiniteElementType::
318  Traits::LocalBasisType::Traits::RangeType RangeType;
319  typedef typename LFSV::Traits::FiniteElementType::
320  Traits::LocalBasisType::Traits::JacobianType JacobianType;
321 
322  // dimensions
323  const int dim = IG::dimension;
324  const int dimw = IG::dimensionworld;
325 
326  // select quadrature rule
327  Dune::GeometryType gtface = ig.geometryInInside().type();
328  const int qorder = std::max ( 2 * ( (int)lfsu.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
329  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
330 
331  // evaluate boundary condition type
332  // Dirichlet boundary condition
333  if( bctype.isDirichlet( ig, rule.begin()->position() ) )
334  {
335  // center in face's reference element
336  const Dune::FieldVector<DF,IG::dimension-1>& face_center =
337  Dune::ReferenceElements<DF,IG::dimension-1>::
338  general(ig.geometry().type()).position(0,0);
339  // outer normal, assuming it is constant over whole intersection
340  const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
341 
342  // evaluate diffusion tensor at cell center, assume it is constant over elements
343  const Dune::FieldVector<DF,IG::dimension>
344  localcenter = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
345  typename K::Traits::RangeType tensor(0.0);
346  k.evaluate(*ig.inside(),localcenter,tensor);
347 
348  // penalty weight for NIPG / SIPG
349  RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
350 
351  // loop over quadrature points and integrate u * phi
352  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
353  {
354  // position of quadrature point in local coordinates of element
355  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
356 
357  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
358  std::vector<JacobianType> js(lfsv.size());
359  lfsv.finiteElement().localBasis().evaluateJacobian(local,js);
360 
361  // transform gradient to real element
362  typename IG::Entity::Geometry::JacobianInverseTransposed jac;
363  jac = ig.inside()->geometry().jacobianInverseTransposed(local);
364  std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsv.size());
365  for (size_t i=0; i<lfsv.size(); i++)
366  {
367  gradphi[i] = 0.0;
368  jac.umv(js[i][0],gradphi[i]);
369  }
370 
371  // evaluate test shape functions
372  std::vector<RangeType> phi(lfsv.size());
373  lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
374 
375  // compute gradient of u
376  Dune::FieldVector<RF,dim> gradu(0.0);
377  for (size_t i=0; i<lfsu.size(); i++)
378  {
379  gradu.axpy(x(lfsu,i),gradphi[i]);
380  }
381 
382  // compute K * gradient of u
383  Dune::FieldVector<RF,dim> Kgradu(0.0);
384  tensor.umv(gradu,Kgradu);
385 
386  // evaluate u
387  RF u=0.0;
388  for (size_t i=0; i<lfsu.size(); i++)
389  {
390  u += x(lfsu,i)*phi[i];
391  }
392 
393  // integrate u
394  RF factor = it->weight()*ig.geometry().integrationElement(it->position());
395  for (size_t i=0; i<lfsv.size(); i++)
396  {
397  // Left hand side of NIPG / SIPG penalty term: sigma/|gamma|^beta*u*v
398  r.accumulate( lfsv, i, penalty_weight*u*phi[i]*factor );
399  // epsilon*K*gradient of v*my*u - v*K*gradient of u*my
400  Dune::FieldVector<RF,dim> Kgradv(0.0);
401  tensor.umv(gradphi[i],Kgradv);
402  r.accumulate( lfsv, i, epsilon * (Kgradv*normal)*u*factor );
403  r.accumulate( lfsv, i, - phi[i]*(Kgradu*normal)*factor );
404  }
405  }
406  }
407  }
408 #endif
409 
410  // volume integral depending only on test functions,
411  // contains f on the right hand side
412  template<typename EG, typename LFSV, typename R>
413  void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
414  {
415  // domain and range field type
416  typedef typename LFSV::Traits::FiniteElementType::
417  Traits::LocalBasisType::Traits::DomainFieldType DF;
418  typedef typename LFSV::Traits::FiniteElementType::
419  Traits::LocalBasisType::Traits::RangeFieldType RF;
420  typedef typename LFSV::Traits::FiniteElementType::
421  Traits::LocalBasisType::Traits::RangeType RangeType;
422 
423  // dimensions
424  const int dim = EG::Geometry::dimension;
425 
426  // select quadrature rule
427  Dune::GeometryType gt = eg.geometry().type();
428  const int qorder = std::max ( 2 * ( (int)lfsv.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
429  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
430 
431  // loop over quadrature points
432  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
433  {
434  // evaluate shape functions
435  std::vector<RangeType> phi(lfsv.size());
436  lfsv.finiteElement().localBasis().evaluateFunction(it->position(),phi);
437 
438  // evaluate right hand side parameter function
439  typename F::Traits::RangeType y;
440  f.evaluate(eg.entity(),it->position(),y);
441 
442  // integrate f
443  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
444  for (size_t i=0; i<lfsv.size(); i++)
445  {
446  // f*v
447  r.accumulate( lfsv, i, - y*phi[i]*factor );
448  }
449  }
450  }
451 
452  // boundary integral independent of ansatz functions,
453  // Neumann and Dirichlet boundary conditions, DG penalty term's right hand side
454  template<typename IG, typename LFSV, typename R>
455  void lambda_boundary (const IG& ig, const LFSV& lfsv, R& r) const
456  {
457  // domain and range field type
458  typedef typename LFSV::Traits::FiniteElementType::
459  Traits::LocalBasisType::Traits::DomainFieldType DF;
460  typedef typename LFSV::Traits::FiniteElementType::
461  Traits::LocalBasisType::Traits::RangeFieldType RF;
462  typedef typename LFSV::Traits::FiniteElementType::
463  Traits::LocalBasisType::Traits::RangeType RangeType;
464  typedef typename LFSV::Traits::FiniteElementType::
465  Traits::LocalBasisType::Traits::JacobianType JacobianType;
466 
467  // dimensions
468  const int dim = IG::dimension;
469  const int dimw = IG::dimensionworld;
470 
471  // select quadrature rule
472  Dune::GeometryType gtface = ig.geometryInInside().type();
473  const int qorder = std::max ( 2 * ( (int)lfsv.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
474  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
475 
476  // evaluate boundary condition type
477  // Neumann boundary condition
478  if( bctype.isNeumann( ig, rule.begin()->position() ) )
479  {
480  // loop over quadrature points and integrate normal flux
481  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
482  {
483  // position of quadrature point in local coordinates of element
484  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
485 
486  // evaluate test shape functions
487  std::vector<RangeType> phi(lfsv.size());
488  lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
489 
490  // evaluate flux boundary condition
491  typename J::Traits::RangeType y(0.0);
492  j.evaluate(*(ig.inside()),local,y);
493 
494  // integrate J
495  RF factor = it->weight()*ig.geometry().integrationElement(it->position());
496  for (size_t i=0; i<lfsv.size(); i++)
497  {
498  // Neumann boundary condition: j*v
499  r.accumulate( lfsv, i, y*phi[i]*factor );
500  }
501  }
502  }
503  // Dirichlet boundary condition
504  else if( bctype.isDirichlet( ig, rule.begin()->position() ) )
505  {
506  /*
507  !!!!!!!! TODO: Warum normale am face center? !!!!!!
508  */
509  // center in face's reference element
510  const Dune::FieldVector<DF,IG::dimension-1>& face_center =
511  Dune::ReferenceElements<DF,IG::dimension-1>::
512  general(ig.geometry().type()).position(0,0);
513  // outer normal, assuming it is constant over whole intersection
514  const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
515  // evaluate diffusion tensor at cell center, assume it is constant over elements
516  const Dune::FieldVector<DF,IG::dimension>
517  localcenter = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
518  typename K::Traits::RangeType tensor(0.0);
519  k.evaluate(*ig.inside(),localcenter,tensor);
520  // penalty weight for NIPG / SIPG
521  RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
522 
523  // loop over quadrature points and integrate g * phi
524  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
525  {
526  // position of quadrature point in local coordinates of element
527  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
528 
529  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
530  std::vector<JacobianType> js(lfsv.size());
531  lfsv.finiteElement().localBasis().evaluateJacobian(local,js);
532 
533  // transform gradient to real element
534  typename IG::Entity::Geometry::JacobianInverseTransposed jac;
535  jac = ig.inside()->geometry().jacobianInverseTransposed(local);
536  std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsv.size());
537  for (size_t i=0; i<lfsv.size(); i++)
538  {
539  gradphi[i] = 0.0;
540  jac.umv(js[i][0],gradphi[i]);
541  }
542 
543  // evaluate test shape functions
544  std::vector<RangeType> phi(lfsv.size());
545  lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
546 
547  // evaluate Dirichlet boundary condition
548  typename G::Traits::RangeType y = 0;
549  g.evaluate(*(ig.inside()),local,y);
550 
551  // integrate G
552  RF factor = it->weight()*ig.geometry().integrationElement(it->position());
553  for (size_t i=0; i<lfsv.size(); i++)
554  {
555  // Right hand side of NIPG / SIPG penalty term: -sigma / |gamma|^beta * g
556  r.accumulate( lfsv, i, -penalty_weight * y * phi[i] * factor );
557  // Dirichlet boundary: -epsilon*K*gradient of v*my*g
558  Dune::FieldVector<RF,dim> Kgradv(0.0);
559  tensor.umv(gradphi[i],Kgradv);
560  r.accumulate( lfsv, i, -epsilon * (Kgradv*normal)*y*factor );
561  }
562  }
563  }
564  else
565  {
566  DUNE_THROW( Dune::Exception, "Unrecognized or unsupported boundary condition type!" );
567  }
568  }
569 
570 #ifndef NumericalJacobianX
571  // jacobian of volume term
572  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
573  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
574  M& mat) const
575  {
576  // domain and range field type
577  typedef typename LFSU::Traits::FiniteElementType::
578  Traits::LocalBasisType::Traits::DomainFieldType DF;
579  typedef typename LFSU::Traits::FiniteElementType::
580  Traits::LocalBasisType::Traits::RangeFieldType RF;
581  typedef typename LFSU::Traits::FiniteElementType::
582  Traits::LocalBasisType::Traits::JacobianType JacobianType;
583 
584  // dimensions
585  const int dim = EG::Geometry::dimension;
586 
587  // select quadrature rule
588  Dune::GeometryType gt = eg.geometry().type();
589  const int qorder = std::max ( 2 * ( (int)lfsu.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
590  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
591 
592  // evaluate diffusion tensor at cell center, assume it is constant over elements
593  typename K::Traits::RangeType tensor;
594  Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
595  k.evaluate(eg.entity(),localcenter,tensor);
596 
597  // loop over quadrature points
598  for (typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
599  {
600  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
601  std::vector<JacobianType> js(lfsu.size());
602  lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
603 
604  // transform gradient to real element
605  typename EG::Geometry::JacobianInverseTransposed jac;
606  jac = eg.geometry().jacobianInverseTransposed(it->position());
607  std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
608  for (typename LFSU::Traits::SizeType i=0; i<lfsu.size(); i++)
609  {
610  gradphi[i] = 0.0;
611  jac.umv(js[i][0],gradphi[i]);
612  }
613 
614  // compute K * gradient of shape functions
615  std::vector<Dune::FieldVector<RF,dim> > Kgradphi(lfsu.size());
616  for (typename LFSU::Traits::SizeType i=0; i<lfsu.size(); i++)
617  {
618  tensor.mv(gradphi[i],Kgradphi[i]);
619  }
620 
621  // integrate (K grad phi_j)*grad phi_i
622  RF factor = it->weight() * eg.geometry().integrationElement(it->position());
623  for (typename LFSU::Traits::SizeType j=0; j<lfsu.size(); j++)
624  {
625  for (typename LFSU::Traits::SizeType i=0; i<lfsu.size(); i++)
626  {
627  // K*gradient of phi_j*K*gradient of phi_i
628  mat.accumulate( lfsu, i, lfsu, j, ( Kgradphi[j]*gradphi[i])*factor );
629  }
630  }
631  }
632  }
633 
634  // jacobian of skeleton term
635  template<typename IG, typename LFSU, typename X, typename LFSV, typename M>
636  void jacobian_skeleton (const IG& ig,
637  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
638  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
639  M& mat_ss, M& mat_sn,
640  M& mat_ns, M& mat_nn) const
641  {
642  // domain and range field type
643  typedef typename LFSU::Traits::FiniteElementType::
644  Traits::LocalBasisType::Traits::DomainFieldType DF;
645  typedef typename LFSU::Traits::FiniteElementType::
646  Traits::LocalBasisType::Traits::RangeFieldType RF;
647  typedef typename LFSV::Traits::FiniteElementType::
648  Traits::LocalBasisType::Traits::RangeType RangeType;
649  typedef typename LFSV::Traits::FiniteElementType::
650  Traits::LocalBasisType::Traits::JacobianType JacobianType;
651 
652  // dimensions
653  const int dim = IG::dimension;
654  const int dimw = IG::dimensionworld;
655 
656  // select quadrature rule
657  Dune::GeometryType gtface = ig.geometryInInside().type();
658  const int qorder = std::max( 0, std::max(
659  2 * ( (int)lfsu_s.finiteElement().localBasis().order() - 1 ),
660  2 * ( (int)lfsu_n.finiteElement().localBasis().order() - 1 ))) + superintegration_order;
661  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
662 
663  // center in face's reference element
664  const Dune::FieldVector<DF,IG::dimension-1>& face_center =
665  Dune::ReferenceElements<DF,IG::dimension-1>::
666  general(ig.geometry().type()).position(0,0);
667  const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
668 
669  // evaluate diffusion tensor at cell center, assume it is constant over elements
670  const Dune::FieldVector<DF,IG::dimension>&
671  inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
672  const Dune::FieldVector<DF,IG::dimension>&
673  outside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.outside()->type()).position(0,0);
674  typename K::Traits::RangeType permeability_s(0.0);
675  typename K::Traits::RangeType permeability_n(0.0);
676  k.evaluate(*(ig.inside()),inside_local,permeability_s);
677  k.evaluate(*(ig.outside()),outside_local,permeability_n);
678 
679  // penalty weight for NIPG / SIPG
680  RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
681 
682  // loop over quadrature points
683  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
684  {
685  // position of quadrature point in local coordinates of element
686  Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(it->position());
687  Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(it->position());
688 
689  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
690  std::vector<JacobianType> js_s(lfsv_s.size());
691  lfsv_s.finiteElement().localBasis().evaluateJacobian(local_s,js_s);
692  std::vector<JacobianType> js_n(lfsv_n.size());
693  lfsv_n.finiteElement().localBasis().evaluateJacobian(local_n,js_n);
694 
695  // transform gradient to real element
696  typename IG::Entity::Geometry::JacobianInverseTransposed jac_s;
697  jac_s = ig.inside()->geometry().jacobianInverseTransposed(local_s);
698  std::vector<Dune::FieldVector<RF,dim> > gradphi_s(lfsv_s.size());
699  for (size_t i=0; i<lfsv_s.size(); i++)
700  {
701  gradphi_s[i] = 0.0;
702  jac_s.umv(js_s[i][0],gradphi_s[i]);
703  }
704  typename IG::Entity::Geometry::JacobianInverseTransposed jac_n;
705  jac_n = ig.outside()->geometry().jacobianInverseTransposed(local_n);
706  std::vector<Dune::FieldVector<RF,dim> > gradphi_n(lfsv_n.size());
707  for (size_t i=0; i<lfsv_n.size(); i++)
708  {
709  gradphi_n[i] = 0.0;
710  jac_n.umv(js_n[i][0],gradphi_n[i]);
711  }
712 
713  // evaluate test shape functions
714  std::vector<RangeType> phi_s(lfsv_s.size());
715  lfsv_s.finiteElement().localBasis().evaluateFunction(local_s,phi_s);
716  std::vector<RangeType> phi_n(lfsv_n.size());
717  lfsv_n.finiteElement().localBasis().evaluateFunction(local_n,phi_n);
718 
719  // compute gradient of u
720  Dune::FieldVector<RF,dim> gradu_s(0.0);
721  for (size_t i=0; i<lfsu_s.size(); i++)
722  {
723  gradu_s.axpy(x_s(lfsu_s,i),gradphi_s[i]);
724  }
725  Dune::FieldVector<RF,dim> gradu_n(0.0);
726  for (size_t i=0; i<lfsu_n.size(); i++)
727  {
728  gradu_n.axpy(x_n(lfsu_n,i),gradphi_n[i]);
729  }
730 
731  // compute K * gradient of u
732  Dune::FieldVector<RF,dim> kgradu_s(0.0);
733  permeability_s.umv(gradu_s,kgradu_s);
734  Dune::FieldVector<RF,dim> kgradu_n(0.0);
735  permeability_n.umv(gradu_n,kgradu_n);
736 
737  // evaluate u in both elements self and neighbor
738  RF u_s = 0.0;
739  for (size_t i=0; i<lfsu_s.size(); i++)
740  {
741  u_s += x_s(lfsu_n,i)*phi_s[i];
742  }
743  RF u_n = 0.0;
744  for (size_t i=0; i<lfsu_n.size(); i++)
745  {
746  u_n += x_n(lfsu_n,i)*phi_n[i];
747  }
748 
749  // average on intersection of K * grad v * normal
750  std::vector<Dune::FieldVector<RF,dim> > kgradphi_s(lfsu_s.size());
751  std::vector<Dune::FieldVector<RF,dim> > kgradphi_n(lfsu_n.size());
752  for (size_t i=0; i<lfsu_s.size(); i++)
753  {
754  permeability_s.mv(gradphi_s[i],kgradphi_s[i]);
755  }
756  for (size_t i=0; i<lfsu_n.size(); i++)
757  {
758  permeability_n.mv(gradphi_n[i],kgradphi_n[i]);
759  }
760 
761  // integrate what needed
762  RF factor = it->weight()*ig.geometry().integrationElement(it->position());
763  for (typename LFSU::Traits::SizeType j=0; j<lfsu_s.size(); j++)
764  {
765  for (typename LFSU::Traits::SizeType i=0; i<lfsu_s.size(); i++)
766  {
767  // epsilon*(K*gradient of phi_s_j + K*gradient of phi_n_j)/2*(phi_s_i - phi_n_i) - (phi_s_j - phi_n_j)*(K*gradient of phi_s_i + K*gradient of phi_n_i)/2
768  mat_ss.accumulate( lfsu_s, i, lfsu_s, j, (epsilon*0.5*(kgradphi_s[i]*normal)*(phi_s[j]) - (phi_s[i])*0.5*(kgradphi_s[j]*normal) ) * factor );
769  // NIPG / SIPG term: (phi_n_j - phi_s_j)*(phi_n_i - phi_s_i)
770  mat_ss.accumulate( lfsu_s, i, lfsu_s, j, penalty_weight*(phi_s[j])*(phi_s[i]) * factor );
771  }
772  for (typename LFSU::Traits::SizeType i=0; i<lfsu_n.size(); i++)
773  {
774  // epsilon*(K*gradient of phi_s_j + K*gradient of phi_n_j)/2*(phi_s_i - phi_n_i) - (phi_s_j - phi_n_j)*(K*gradient of phi_s_i + K*gradient of phi_n_i)/2
775  mat_ns.accumulate( lfsu_n, i, lfsu_s, j, (epsilon*0.5*(kgradphi_n[i]*normal)*(phi_s[j]) - (-phi_n[i])*0.5*(kgradphi_s[j]*normal) )*factor );
776  // NIPG / SIPG term: (phi_n_j - phi_s_j)*(phi_n_i - phi_s_i)
777  mat_ns.accumulate( lfsu_n, i, lfsu_s, j, penalty_weight*(phi_s[j])*(-phi_n[i]) *factor );
778  }
779  }
780  for (typename LFSU::Traits::SizeType j=0; j<lfsu_n.size(); j++)
781  {
782  for (typename LFSU::Traits::SizeType i=0; i<lfsu_s.size(); i++)
783  {
784  // epsilon*(K*gradient of phi_s_j + K*gradient of phi_n_j)/2*(phi_s_i - phi_n_i) - (phi_s_j - phi_n_j)*(K*gradient of phi_s_i + K*gradient of phi_n_i)/2
785  mat_sn.accumulate( lfsu_s, i, lfsu_n, j, (epsilon*0.5*(kgradphi_s[i]*normal)*(-phi_n[j]) - (phi_s[i])*0.5*(kgradphi_n[j]*normal) )*factor );
786  // NIPG / SIPG term: (phi_n_j - phi_s_j)*(phi_n_i - phi_s_i)
787  mat_sn.accumulate( lfsu_s, i, lfsu_n, j, penalty_weight*(-phi_n[j])*(phi_s[i]) *factor );
788  }
789  for (typename LFSU::Traits::SizeType i=0; i<lfsu_n.size(); i++)
790  {
791  // epsilon*(K*gradient of phi_s_j + K*gradient of phi_n_j)/2*(phi_s_i - phi_n_i) - (phi_s_j - phi_n_j)*(K*gradient of phi_s_i + K*gradient of phi_n_i)/2
792  mat_nn.accumulate( lfsu_s, i, lfsu_n, j, (epsilon*0.5*(kgradphi_n[i]*normal)*(-phi_n[j]) - (-phi_n[i])*0.5*(kgradphi_n[j]*normal) )*factor );
793  // NIPG / SIPG term: (phi_n_j - phi_s_j)*(phi_n_i - phi_s_i)
794  mat_nn.accumulate( lfsu_s, i, lfsu_n, j, penalty_weight*(-phi_n[j])*(-phi_n[i]) *factor );
795  }
796  }
797  }
798  }
799 
800  // jacobian of volume term
801  template<typename IG, typename LFSU, typename X, typename LFSV, typename M>
802  void jacobian_boundary (const IG& ig,
803  const LFSU& lfsu, const X& x, const LFSV& lfsv,
804  M& mat) const
805  {
806  // domain and range field type
807  typedef typename LFSU::Traits::FiniteElementType::
808  Traits::LocalBasisType::Traits::DomainFieldType DF;
809  typedef typename LFSU::Traits::FiniteElementType::
810  Traits::LocalBasisType::Traits::RangeFieldType RF;
811  typedef typename LFSV::Traits::FiniteElementType::
812  Traits::LocalBasisType::Traits::RangeType RangeType;
813  typedef typename LFSV::Traits::FiniteElementType::
814  Traits::LocalBasisType::Traits::JacobianType JacobianType;
815 
816  // dimensions
817  const int dim = IG::dimension;
818  const int dimw = IG::dimensionworld;
819 
820  // select quadrature rule
821  Dune::GeometryType gtface = ig.geometryInInside().type();
822  const int qorder = std::max ( 2 * ( (int)lfsu.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
823  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
824 
825  // evaluate boundary condition type
826  // Dirichlet boundary condition
827  if( bctype.isDirichlet( ig, rule.begin()->position() ) )
828  {
829  // center in face's reference element
830  const Dune::FieldVector<DF,IG::dimension-1>& face_center =
831  Dune::ReferenceElements<DF,IG::dimension-1>::
832  general(ig.geometry().type()).position(0,0);
833  // outer normal, assuming it is constant over whole intersection
834  const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
835 
836  // evaluate diffusion tensor at cell center, assume it is constant over elements
837  const Dune::FieldVector<DF,IG::dimension>
838  localcenter = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
839  typename K::Traits::RangeType tensor(0.0);
840  k.evaluate(*ig.inside(),localcenter,tensor);
841 
842  // penalty weight for NIPG / SIPG
843  RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
844 
845  // loop over quadrature points and integrate u * phi
846  for (typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
847  {
848  // position of quadrature point in local coordinates of element
849  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
850 
851  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
852  std::vector<JacobianType> js(lfsv.size());
853  lfsv.finiteElement().localBasis().evaluateJacobian(local,js);
854 
855  // transform gradient to real element
856  typename IG::Entity::Geometry::JacobianInverseTransposed jac;
857  jac = ig.inside()->geometry().jacobianInverseTransposed(local);
858  std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsv.size());
859  for (size_t i=0; i<lfsv.size(); i++)
860  {
861  gradphi[i] = 0.0;
862  jac.umv(js[i][0],gradphi[i]);
863  }
864 
865  // evaluate test shape functions
866  std::vector<RangeType> phi(lfsv.size());
867  lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
868 
869  // compute gradient of u
870  Dune::FieldVector<RF,dim> gradu(0.0);
871  for (size_t i=0; i<lfsu.size(); i++)
872  {
873  gradu.axpy(x(lfsu,i),gradphi[i]);
874  }
875 
876  // compute K * gradient of v
877  std::vector<Dune::FieldVector<RF,dim> > kgradphi(lfsu.size());
878  for (size_t i=0; i<lfsu.size(); i++)
879  {
880  tensor.mv(gradphi[i],kgradphi[i]);
881  }
882 
883  // integrate
884  RF factor = it->weight()*ig.geometry().integrationElement(it->position());
885  for (typename LFSU::Traits::SizeType j=0; j<lfsu.size(); j++)
886  {
887  for (typename LFSU::Traits::SizeType i=0; i<lfsu.size(); i++)
888  {
889  // epsilon*K*gradient of phi_i*my*phi_j - phi_i*K*gradient of phi_j*my
890  mat.accumulate( lfsu, i, lfsu, j, (epsilon*(kgradphi[i]*normal)*phi[j] - phi[i]*(kgradphi[j]*normal))*factor );
891  // NIPG / SIPG penalty term: sigma / |gamma|^beta *phi_j*phi_i
892  mat.accumulate( lfsu, i, lfsu, j, (penalty_weight*phi[j]*phi[i])*factor );
893  }
894  }
895  }
896  }
897  }
898 #endif
899 
900  private:
901  const K& k;
902  const F& f;
903  const B& bctype;
904  const G& g;
905  const J& j;
906  // values for NIPG / NIPG
907  double epsilon;
908  double sigma;
909  double beta;
910  int superintegration_order; // Quadrature order
911  };
912 
914  } // namespace PDELab
915 } // namespace Dune
916 
917 #endif
Definition: diffusiondg.hh:59
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusiondg.hh:96
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: diffusiondg.hh:455
Definition: diffusiondg.hh:60
sparsity pattern generator
Definition: pattern.hh:13
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: diffusiondg.hh:636
const IG & ig
Definition: constraints.hh:148
static const int dim
Definition: adaptivity.hh:83
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: diffusiondg.hh:573
void jacobian_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: diffusiondg.hh:802
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: diffusiondg.hh:157
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:250
Definition: adaptivity.hh:27
Implement alpha_boundary() based on jacobian_boundary()
Definition: numericalresidual.hh:128
sparsity pattern generator
Definition: pattern.hh:29
Default flags for all local operators.
Definition: flags.hh:18
Definition: diffusiondg.hh:54
Definition: diffusiondg.hh:61
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:31
Definition: diffusiondg.hh:63
Implement alpha_volume() based on jacobian_volume()
Definition: numericalresidual.hh:35
void alpha_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusiondg.hh:310
Definition: diffusiondg.hh:58
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: numericaljacobian.hh:156
Definition: diffusiondg.hh:62
Implement alpha_skeleton() based on jacobian_skeleton()
Definition: numericalresidual.hh:74
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: diffusiondg.hh:413
Definition: diffusiondg.hh:33