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