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