dune-pdelab  2.4-dev
convectiondiffusionfem.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_CONVECTIONDIFFUSIONFEM_HH
3 #define DUNE_PDELAB_CONVECTIONDIFFUSIONFEM_HH
4 
5 #include<vector>
6 
7 #include<dune/common/exceptions.hh>
8 #include<dune/common/fvector.hh>
9 #include<dune/geometry/type.hh>
10 #include<dune/geometry/referenceelements.hh>
11 #include<dune/geometry/quadraturerules.hh>
12 
18 
20 
21 namespace Dune {
22  namespace PDELab {
23 
38  template<typename T, typename FiniteElementMap>
40  public Dune::PDELab::NumericalJacobianApplyVolume<ConvectionDiffusionFEM<T,FiniteElementMap> >,
41  public Dune::PDELab::NumericalJacobianApplyBoundary<ConvectionDiffusionFEM<T,FiniteElementMap> >,
42  public Dune::PDELab::NumericalJacobianVolume<ConvectionDiffusionFEM<T,FiniteElementMap> >,
43  public Dune::PDELab::NumericalJacobianBoundary<ConvectionDiffusionFEM<T,FiniteElementMap> >,
46  public Dune::PDELab::InstationaryLocalOperatorDefaultMethods<typename T::Traits::RangeFieldType>
47  {
48  public:
49  // pattern assembly flags
50  enum { doPatternVolume = true };
51 
52  // residual assembly flags
53  enum { doAlphaVolume = true };
54  enum { doAlphaBoundary = true };
55 
56  ConvectionDiffusionFEM (T& param_, int intorderadd_=0)
57  : param(param_), intorderadd(intorderadd_)
58  {
59  }
60 
61  // volume integral depending on test and ansatz functions
62  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
63  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
64  {
65  // domain and range field type
66  typedef typename LFSU::Traits::FiniteElementType::
67  Traits::LocalBasisType::Traits::DomainFieldType DF;
68  typedef typename LFSU::Traits::FiniteElementType::
69  Traits::LocalBasisType::Traits::RangeFieldType RF;
70  typedef typename LFSU::Traits::FiniteElementType::
71  Traits::LocalBasisType::Traits::JacobianType JacobianType;
72  typedef typename LFSU::Traits::FiniteElementType::
73  Traits::LocalBasisType::Traits::RangeType RangeType;
74 
75  typedef typename LFSU::Traits::SizeType size_type;
76 
77  // dimensions
78  const int dim = EG::Entity::dimension;
79 
80  // select quadrature rule
81  Dune::GeometryType gt = eg.geometry().type();
82  const int intorder = intorderadd+2*lfsu.finiteElement().localBasis().order();
83  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
84 
85  // evaluate diffusion tensor at cell center, assume it is constant over elements
86  typename T::Traits::PermTensorType tensor;
87  Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
88  tensor = param.A(eg.entity(),localcenter);
89 
90  // loop over quadrature points
91  for (const auto& ip : rule)
92  {
93  // evaluate basis functions
94  // std::vector<RangeType> phi(lfsu.size());
95  // lfsu.finiteElement().localBasis().evaluateFunction(ip.position(),phi);
96  const std::vector<RangeType>& phi = cache.evaluateFunction(ip.position(),lfsu.finiteElement().localBasis());
97 
98  // evaluate u
99  RF u=0.0;
100  for (size_type i=0; i<lfsu.size(); i++)
101  u += x(lfsu,i)*phi[i];
102 
103  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
104  // std::vector<JacobianType> js(lfsu.size());
105  // lfsu.finiteElement().localBasis().evaluateJacobian(ip.position(),js);
106  const std::vector<JacobianType>& js = cache.evaluateJacobian(ip.position(),lfsu.finiteElement().localBasis());
107 
108  // transform gradients of shape functions to real element
109  const typename EG::Geometry::JacobianInverseTransposed jac =
110  eg.geometry().jacobianInverseTransposed(ip.position());
111  std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
112  for (size_type i=0; i<lfsu.size(); i++)
113  jac.mv(js[i][0],gradphi[i]);
114 
115  // compute gradient of u
116  Dune::FieldVector<RF,dim> gradu(0.0);
117  for (size_type i=0; i<lfsu.size(); i++)
118  gradu.axpy(x(lfsu,i),gradphi[i]);
119 
120  // compute A * gradient of u
121  Dune::FieldVector<RF,dim> Agradu(0.0);
122  tensor.umv(gradu,Agradu);
123 
124  // evaluate velocity field, sink term and source term
125  typename T::Traits::RangeType b = param.b(eg.entity(),ip.position());
126  typename T::Traits::RangeFieldType c = param.c(eg.entity(),ip.position());
127  typename T::Traits::RangeFieldType f = param.f(eg.entity(),ip.position());
128 
129  // integrate (A grad u)*grad phi_i - u b*grad phi_i + c*u*phi_i
130  RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
131  for (size_type i=0; i<lfsu.size(); i++)
132  r.accumulate(lfsu,i,( Agradu*gradphi[i] - u*(b*gradphi[i]) + (c*u-f)*phi[i] )*factor);
133  }
134  }
135 
136  // jacobian of volume term
137  template<typename EG, typename LFSU, typename X, typename LFSV, typename M>
138  void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
139  M& mat) const
140  {
141  // domain and range field type
142  typedef typename LFSU::Traits::FiniteElementType::
143  Traits::LocalBasisType::Traits::DomainFieldType DF;
144  typedef typename LFSU::Traits::FiniteElementType::
145  Traits::LocalBasisType::Traits::RangeFieldType RF;
146  typedef typename LFSU::Traits::FiniteElementType::
147  Traits::LocalBasisType::Traits::JacobianType JacobianType;
148  typedef typename LFSU::Traits::FiniteElementType::
149  Traits::LocalBasisType::Traits::RangeType RangeType;
150  typedef typename LFSU::Traits::SizeType size_type;
151 
152  // dimensions
153  const int dim = EG::Entity::dimension;
154 
155  // select quadrature rule
156  Dune::GeometryType gt = eg.geometry().type();
157  const int intorder = intorderadd+2*lfsu.finiteElement().localBasis().order();
158  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
159 
160  // evaluate diffusion tensor at cell center, assume it is constant over elements
161  typename T::Traits::PermTensorType tensor;
162  Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
163  tensor = param.A(eg.entity(),localcenter);
164 
165  // loop over quadrature points
166  for (const auto& ip : rule)
167  {
168  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
169  // std::vector<JacobianType> js(lfsu.size());
170  // lfsu.finiteElement().localBasis().evaluateJacobian(ip.position(),js);
171  const std::vector<JacobianType>& js = cache.evaluateJacobian(ip.position(),lfsu.finiteElement().localBasis());
172 
173  // transform gradient to real element
174  const typename EG::Geometry::JacobianInverseTransposed jac
175  = eg.geometry().jacobianInverseTransposed(ip.position());
176  std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
177  std::vector<Dune::FieldVector<RF,dim> > Agradphi(lfsu.size());
178  for (size_type i=0; i<lfsu.size(); i++)
179  {
180  jac.mv(js[i][0],gradphi[i]);
181  tensor.mv(gradphi[i],Agradphi[i]);
182  }
183 
184  // evaluate basis functions
185  // std::vector<RangeType> phi(lfsu.size());
186  // lfsu.finiteElement().localBasis().evaluateFunction(ip.position(),phi);
187  const std::vector<RangeType>& phi = cache.evaluateFunction(ip.position(),lfsu.finiteElement().localBasis());
188 
189  // evaluate velocity field, sink term and source te
190  typename T::Traits::RangeType b = param.b(eg.entity(),ip.position());
191  typename T::Traits::RangeFieldType c = param.c(eg.entity(),ip.position());
192 
193  // integrate (A grad phi_j)*grad phi_i - phi_j b*grad phi_i + c*phi_j*phi_i
194  RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
195  for (size_type j=0; j<lfsu.size(); j++)
196  for (size_type i=0; i<lfsu.size(); i++)
197  mat.accumulate(lfsu,i,lfsu,j,( Agradphi[j]*gradphi[i]-phi[j]*(b*gradphi[i])+c*phi[j]*phi[i] )*factor);
198  }
199  }
200 
201  // boundary integral
202  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
203  void alpha_boundary (const IG& ig,
204  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
205  R& r_s) const
206  {
207  // domain and range field type
208  typedef typename LFSV::Traits::FiniteElementType::
209  Traits::LocalBasisType::Traits::DomainFieldType DF;
210  typedef typename LFSV::Traits::FiniteElementType::
211  Traits::LocalBasisType::Traits::RangeFieldType RF;
212  typedef typename LFSV::Traits::FiniteElementType::
213  Traits::LocalBasisType::Traits::RangeType RangeType;
214 
215  typedef typename LFSV::Traits::SizeType size_type;
216 
217  // dimensions
218  const int dim = IG::dimension;
219 
220  // get cell entity
221  auto inside_cell = ig.inside();
222 
223  // evaluate boundary condition type
224  Dune::GeometryType gtface = ig.geometryInInside().type();
225  Dune::FieldVector<DF,dim-1> facecenterlocal = Dune::ReferenceElements<DF,dim-1>::general(gtface).position(0,0);
227  bctype = param.bctype(ig.intersection(),facecenterlocal);
228 
229  // skip rest if we are on Dirichlet boundary
231 
232  // select quadrature rule
233  const int intorder = intorderadd+2*lfsu_s.finiteElement().localBasis().order();
234  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
235 
236  // loop over quadrature points and integrate normal flux
237  for (const auto& ip : rule)
238  {
239  // position of quadrature point in local coordinates of element
240  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
241 
242  // evaluate shape functions (assume Galerkin method)
243  // std::vector<RangeType> phi(lfsu_s.size());
244  // lfsu_s.finiteElement().localBasis().evaluateFunction(local,phi);
245  const std::vector<RangeType>& phi = cache.evaluateFunction(local,lfsu_s.finiteElement().localBasis());
246 
248  {
249  // evaluate flux boundary condition
250  typename T::Traits::RangeFieldType j = param.j(ig.intersection(),ip.position());
251 
252  // integrate j
253  RF factor = ip.weight()*ig.geometry().integrationElement(ip.position());
254  for (size_type i=0; i<lfsu_s.size(); i++)
255  r_s.accumulate(lfsu_s,i,j*phi[i]*factor);
256  }
257 
259  {
260  // evaluate u
261  RF u=0.0;
262  for (size_type i=0; i<lfsu_s.size(); i++)
263  u += x_s(lfsu_s,i)*phi[i];
264 
265  // evaluate velocity field and outer unit normal
266  typename T::Traits::RangeType b = param.b(inside_cell,local);
267  const Dune::FieldVector<DF,dim> n = ig.unitOuterNormal(ip.position());
268 
269  // evaluate outflow boundary condition
270  typename T::Traits::RangeFieldType o = param.o(ig.intersection(),ip.position());
271 
272  // integrate o
273  RF factor = ip.weight()*ig.geometry().integrationElement(ip.position());
274  for (size_type i=0; i<lfsu_s.size(); i++)
275  r_s.accumulate(lfsu_s,i,( (b*n)*u + o)*phi[i]*factor);
276  }
277  }
278  }
279 
280  // jacobian contribution from boundary
281  template<typename IG, typename LFSU, typename X, typename LFSV, typename M>
282  void jacobian_boundary (const IG& ig,
283  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
284  M& mat_s) const
285  {
286  // domain and range field type
287  typedef typename LFSV::Traits::FiniteElementType::
288  Traits::LocalBasisType::Traits::DomainFieldType DF;
289  typedef typename LFSV::Traits::FiniteElementType::
290  Traits::LocalBasisType::Traits::RangeFieldType RF;
291  typedef typename LFSV::Traits::FiniteElementType::
292  Traits::LocalBasisType::Traits::RangeType RangeType;
293 
294  typedef typename LFSV::Traits::SizeType size_type;
295 
296  // dimensions
297  const int dim = IG::dimension;
298 
299  // get cell entity
300  auto inside_cell = ig.inside();
301 
302  // evaluate boundary condition type
303  Dune::GeometryType gtface = ig.geometryInInside().type();
304  Dune::FieldVector<DF,dim-1> facecenterlocal = Dune::ReferenceElements<DF,dim-1>::general(gtface).position(0,0);
306  bctype = param.bctype(ig.intersection(),facecenterlocal);
307 
308  // skip rest if we are on Dirichlet boundary
311 
312  // select quadrature rule
313  const int intorder = intorderadd+2*lfsu_s.finiteElement().localBasis().order();
314  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
315 
316  // loop over quadrature points and integrate normal flux
317  for (const auto& ip : rule)
318  {
319  // position of quadrature point in local coordinates of element
320  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
321 
322  // evaluate shape functions (assume Galerkin method)
323  // std::vector<RangeType> phi(lfsu_s.size());
324  // lfsu_s.finiteElement().localBasis().evaluateFunction(local,phi);
325  const std::vector<RangeType>& phi = cache.evaluateFunction(local,lfsu_s.finiteElement().localBasis());
326 
327  // evaluate velocity field and outer unit normal
328  typename T::Traits::RangeType b = param.b(inside_cell,local);
329  const Dune::FieldVector<DF,dim> n = ig.unitOuterNormal(ip.position());
330 
331  // integrate
332  RF factor = ip.weight()*ig.geometry().integrationElement(ip.position());
333  for (size_type j=0; j<lfsu_s.size(); j++)
334  for (size_type i=0; i<lfsu_s.size(); i++)
335  mat_s.accumulate(lfsu_s,i,lfsu_s,j,(b*n)*phi[j]*phi[i]*factor);
336  }
337  }
338 
339 
341  void setTime (double t)
342  {
343  param.setTime(t);
344  }
345 
346  private:
347  T& param;
348  int intorderadd;
349  typedef typename FiniteElementMap::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
351  };
352 
353 
354 
370  template<typename T>
373  {
374  enum { dim = T::Traits::GridViewType::dimension };
375 
376  typedef typename T::Traits::RangeFieldType Real;
378 
379  public:
380  // pattern assembly flags
381  enum { doPatternVolume = false };
382  enum { doPatternSkeleton = false };
383 
384  // residual assembly flags
385  enum { doAlphaVolume = true };
386  enum { doAlphaSkeleton = true };
387  enum { doAlphaBoundary = true };
388 
391  : param(param_)
392  {}
393 
394  // volume integral depending on test and ansatz functions
395  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
396  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
397  {
398  // domain and range field type
399  typedef typename LFSU::Traits::FiniteElementType::
400  Traits::LocalBasisType::Traits::DomainFieldType DF;
401  typedef typename LFSU::Traits::FiniteElementType::
402  Traits::LocalBasisType::Traits::RangeFieldType RF;
403  typedef typename LFSU::Traits::FiniteElementType::
404  Traits::LocalBasisType::Traits::RangeType RangeType;
405  typedef typename LFSU::Traits::SizeType size_type;
406 
407  // dimensions
408  const int dim = EG::Geometry::dimension;
409  const int intorder = 2*lfsu.finiteElement().localBasis().order();
410 
411  // select quadrature rule
412  Dune::GeometryType gt = eg.geometry().type();
413  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
414 
415  // loop over quadrature points
416  RF sum(0.0);
417  for (const auto& ip : rule)
418  {
419  // evaluate basis functions
420  std::vector<RangeType> phi(lfsu.size());
421  lfsu.finiteElement().localBasis().evaluateFunction(ip.position(),phi);
422 
423  // evaluate u
424  RF u=0.0;
425  for (size_type i=0; i<lfsu.size(); i++)
426  u += x(lfsu,i)*phi[i];
427 
428  // evaluate reaction term
429  typename T::Traits::RangeFieldType c = param.c(eg.entity(),ip.position());
430 
431  // evaluate right hand side parameter function
432  typename T::Traits::RangeFieldType f = param.f(eg.entity(),ip.position());
433 
434  // integrate f^2
435  RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
436  sum += (f*f-c*c*u*u)*factor;
437  }
438 
439  // accumulate cell indicator
440  DF h_T = diameter(eg.geometry());
441  r.accumulate(lfsv,0,h_T*h_T*sum);
442  }
443 
444 
445  // skeleton integral depending on test and ansatz functions
446  // each face is only visited ONCE!
447  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
448  void alpha_skeleton (const IG& ig,
449  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
450  const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
451  R& r_s, R& r_n) const
452  {
453  // domain and range field type
454  typedef typename LFSU::Traits::FiniteElementType::
455  Traits::LocalBasisType::Traits::DomainFieldType DF;
456  typedef typename LFSU::Traits::FiniteElementType::
457  Traits::LocalBasisType::Traits::RangeFieldType RF;
458  typedef typename LFSU::Traits::FiniteElementType::
459  Traits::LocalBasisType::Traits::JacobianType JacobianType;
460  typedef typename LFSU::Traits::SizeType size_type;
461 
462  // dimensions
463  const int dim = IG::dimension;
464 
465  // get cell entities from both sides of the intersection
466  auto inside_cell = ig.inside();
467  auto outside_cell = ig.outside();
468 
469  // evaluate permeability tensors
470  const Dune::FieldVector<DF,dim>&
471  inside_local = Dune::ReferenceElements<DF,dim>::general(inside_cell.type()).position(0,0);
472  const Dune::FieldVector<DF,dim>&
473  outside_local = Dune::ReferenceElements<DF,dim>::general(outside_cell.type()).position(0,0);
474  typename T::Traits::PermTensorType A_s, A_n;
475  A_s = param.A(inside_cell,inside_local);
476  A_n = param.A(outside_cell,outside_local);
477 
478  // select quadrature rule
479  const int intorder = 2*lfsu_s.finiteElement().localBasis().order();
480  Dune::GeometryType gtface = ig.geometryInInside().type();
481  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
482 
483  // transformation
484  typename IG::Entity::Geometry::JacobianInverseTransposed jac;
485 
486  // tensor times normal
487  const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
488  Dune::FieldVector<RF,dim> An_F_s;
489  A_s.mv(n_F,An_F_s);
490  Dune::FieldVector<RF,dim> An_F_n;
491  A_n.mv(n_F,An_F_n);
492 
493  // loop over quadrature points and integrate normal flux
494  RF sum(0.0);
495  for (const auto& ip : rule)
496  {
497  // position of quadrature point in local coordinates of elements
498  Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(ip.position());
499  Dune::FieldVector<DF,dim> iplocal_n = ig.geometryInOutside().global(ip.position());
500 
501  // evaluate gradient of basis functions
502  std::vector<JacobianType> gradphi_s(lfsu_s.size());
503  lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
504  std::vector<JacobianType> gradphi_n(lfsu_n.size());
505  lfsu_n.finiteElement().localBasis().evaluateJacobian(iplocal_n,gradphi_n);
506 
507  // transform gradients of shape functions to real element
508  jac = inside_cell.geometry().jacobianInverseTransposed(iplocal_s);
509  std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
510  for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
511  jac = outside_cell.geometry().jacobianInverseTransposed(iplocal_n);
512  std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
513  for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
514 
515  // compute gradient of u
516  Dune::FieldVector<RF,dim> gradu_s(0.0);
517  for (size_type i=0; i<lfsu_s.size(); i++)
518  gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
519  Dune::FieldVector<RF,dim> gradu_n(0.0);
520  for (size_type i=0; i<lfsu_n.size(); i++)
521  gradu_n.axpy(x_n(lfsu_n,i),tgradphi_n[i]);
522 
523  // integrate
524  RF factor = ip.weight() * ig.geometry().integrationElement(ip.position());
525  RF jump = (An_F_s*gradu_s)-(An_F_n*gradu_n);
526  sum += 0.25*jump*jump*factor;
527  }
528 
529  // accumulate indicator
530  // DF h_T = diameter(ig.geometry());
531  DF h_T = std::max(diameter(inside_cell.geometry()),diameter(outside_cell.geometry()));
532  r_s.accumulate(lfsv_s,0,h_T*sum);
533  r_n.accumulate(lfsv_n,0,h_T*sum);
534  }
535 
536 
537  // boundary integral depending on test and ansatz functions
538  // We put the Dirchlet evaluation also in the alpha term to save some geometry evaluations
539  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
540  void alpha_boundary (const IG& ig,
541  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
542  R& r_s) const
543  {
544  // domain and range field type
545  typedef typename LFSU::Traits::FiniteElementType::
546  Traits::LocalBasisType::Traits::DomainFieldType DF;
547  typedef typename LFSU::Traits::FiniteElementType::
548  Traits::LocalBasisType::Traits::RangeFieldType RF;
549  typedef typename LFSU::Traits::FiniteElementType::
550  Traits::LocalBasisType::Traits::JacobianType JacobianType;
551  typedef typename LFSU::Traits::SizeType size_type;
552 
553  // dimensions
554  const int dim = IG::dimension;
555 
556  auto inside_cell = ig.inside();
557  // evaluate permeability tensors
558  const Dune::FieldVector<DF,dim>&
559  inside_local = Dune::ReferenceElements<DF,dim>::general(inside_cell.type()).position(0,0);
560  typename T::Traits::PermTensorType A_s;
561  A_s = param.A(inside_cell,inside_local);
562  const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
563  Dune::FieldVector<RF,dim> An_F_s;
564  A_s.mv(n_F,An_F_s);
565 
566  // select quadrature rule
567  const int intorder = 2*lfsu_s.finiteElement().localBasis().order();
568  Dune::GeometryType gtface = ig.geometryInInside().type();
569  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
570 
571  // transformation
572  typename IG::Entity::Geometry::JacobianInverseTransposed jac;
573 
574  // evaluate boundary condition
575  const Dune::FieldVector<DF,dim-1>
576  face_local = Dune::ReferenceElements<DF,dim-1>::general(gtface).position(0,0);
577  BCType bctype = param.bctype(ig.intersection(),face_local);
579  return;
580 
581  // loop over quadrature points and integrate normal flux
582  RF sum(0.0);
583  for (const auto& ip : rule)
584  {
585  // position of quadrature point in local coordinates of elements
586  Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(ip.position());
587 
588  // evaluate gradient of basis functions
589  std::vector<JacobianType> gradphi_s(lfsu_s.size());
590  lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
591 
592  // transform gradients of shape functions to real element
593  jac = inside_cell.geometry().jacobianInverseTransposed(iplocal_s);
594  std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
595  for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
596 
597  // compute gradient of u
598  Dune::FieldVector<RF,dim> gradu_s(0.0);
599  for (size_type i=0; i<lfsu_s.size(); i++)
600  gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
601 
602  // evaluate flux boundary condition
603  RF j = param.j(ig.intersection(),ip.position());
604 
605  // integrate
606  RF factor = ip.weight() * ig.geometry().integrationElement(ip.position());
607  RF jump = j+(An_F_s*gradu_s);
608  sum += jump*jump*factor;
609  }
610 
611  // accumulate indicator
612  //DF h_T = diameter(ig.geometry());
613  DF h_T = diameter(inside_cell.geometry());
614  r_s.accumulate(lfsv_s,0,h_T*sum);
615  }
616 
617  private:
618  T& param; // two phase parameter class
619 
620  template<class GEO>
621  typename GEO::ctype diameter (const GEO& geo) const
622  {
623  typedef typename GEO::ctype DF;
624  DF hmax = -1.0E00;
625  const int dim = GEO::coorddimension;
626  for (int i=0; i<geo.corners(); i++)
627  {
628  Dune::FieldVector<DF,dim> xi = geo.corner(i);
629  for (int j=i+1; j<geo.corners(); j++)
630  {
631  Dune::FieldVector<DF,dim> xj = geo.corner(j);
632  xj -= xi;
633  hmax = std::max(hmax,xj.two_norm());
634  }
635  }
636  return hmax;
637  }
638 
639  };
640 
641 
657  template<typename T>
661  {
662  enum { dim = T::Traits::GridViewType::dimension };
663 
664  typedef typename T::Traits::RangeFieldType Real;
665  typedef typename ConvectionDiffusionBoundaryConditions::Type BCType;
666 
667  public:
668  // pattern assembly flags
669  enum { doPatternVolume = false };
670  enum { doPatternSkeleton = false };
671 
672  // residual assembly flags
673  enum { doAlphaVolume = true };
674 
676  // supply time step from implicit Euler scheme
677  ConvectionDiffusionTemporalResidualEstimator1 (T& param_, double time_, double dt_)
678  : param(param_), time(time_), dt(dt_), cmax(0)
679  {}
680 
681  // volume integral depending on test and ansatz functions
682  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
683  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
684  {
685  // domain and range field type
686  typedef typename LFSU::Traits::FiniteElementType::
687  Traits::LocalBasisType::Traits::DomainFieldType DF;
688  typedef typename LFSU::Traits::FiniteElementType::
689  Traits::LocalBasisType::Traits::RangeFieldType RF;
690  typedef typename LFSU::Traits::FiniteElementType::
691  Traits::LocalBasisType::Traits::RangeType RangeType;
692  typedef typename LFSU::Traits::SizeType size_type;
693 
694  // dimensions
695  const int dim = EG::Geometry::dimension;
696  const int intorder = 2*lfsu.finiteElement().localBasis().order();
697 
698  // select quadrature rule
699  Dune::GeometryType gt = eg.geometry().type();
700  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
701 
702  // loop over quadrature points
703  RF sum(0.0);
704  RF fsum_up(0.0);
705  RF fsum_mid(0.0);
706  RF fsum_down(0.0);
707  for (const auto& ip : rule)
708  {
709  // evaluate basis functions
710  std::vector<RangeType> phi(lfsu.size());
711  lfsu.finiteElement().localBasis().evaluateFunction(ip.position(),phi);
712 
713  // evaluate u
714  RF u=0.0;
715  for (size_type i=0; i<lfsu.size(); i++)
716  u += x(lfsu,i)*phi[i];
717 
718  // integrate f^2
719  RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
720  sum += u*u*factor;
721 
722  // evaluate right hand side parameter function
723  param.setTime(time);
724  typename T::Traits::RangeFieldType f_down = param.f(eg.entity(),ip.position());
725  param.setTime(time+0.5*dt);
726  typename T::Traits::RangeFieldType f_mid = param.f(eg.entity(),ip.position());
727  param.setTime(time+dt);
728  typename T::Traits::RangeFieldType f_up = param.f(eg.entity(),ip.position());
729  RF f_average = (1.0/6.0)*f_down + (2.0/3.0)*f_mid + (1.0/6.0)*f_up;
730 
731  // integrate f-f_average
732  fsum_down += (f_down-f_average)*(f_down-f_average)*factor;
733  fsum_mid += (f_mid-f_average)*(f_mid-f_average)*factor;
734  fsum_up += (f_up-f_average)*(f_up-f_average)*factor;
735  }
736 
737  // accumulate cell indicator
738  DF h_T = diameter(eg.geometry());
739  r.accumulate(lfsv,0,(h_T*h_T/dt)*sum); // h^2*k_n||jump/k_n||^2
740  r.accumulate(lfsv,0,h_T*h_T * dt*((1.0/6.0)*fsum_down+(2.0/3.0)*fsum_mid+(1.0/6.0)*fsum_up) ); // h^2*||f-time_average(f)||^2_0_s_t
741  }
742 
743  void clearCmax ()
744  {
745  cmax = 0;
746  }
747 
748  double getCmax () const
749  {
750  return cmax;
751  }
752 
753  private:
754  T& param; // two phase parameter class
755  double time;
756  double dt;
757  mutable double cmax;
758 
759  template<class GEO>
760  typename GEO::ctype diameter (const GEO& geo) const
761  {
762  typedef typename GEO::ctype DF;
763  DF hmax = -1.0E00;
764  const int dim = GEO::coorddimension;
765  for (int i=0; i<geo.corners(); i++)
766  {
767  Dune::FieldVector<DF,dim> xi = geo.corner(i);
768  for (int j=i+1; j<geo.corners(); j++)
769  {
770  Dune::FieldVector<DF,dim> xj = geo.corner(j);
771  xj -= xi;
772  hmax = std::max(hmax,xj.two_norm());
773  }
774  }
775  return hmax;
776  }
777 
778  };
779 
780  // a functor that can be used to evaluate rhs parameter function in interpolate
781  template<typename T, typename EG>
783  {
784  public:
785  CD_RHS_LocalAdapter (const T& t_, const EG& eg_) : t(t_), eg(eg_)
786  {}
787 
788  template<typename X, typename Y>
789  inline void evaluate (const X& x, Y& y) const
790  {
791  y[0] = t.f(eg.entity(),x);
792  }
793 
794  private:
795  const T& t;
796  const EG& eg;
797  };
798 
814  template<typename T>
818  {
819  enum { dim = T::Traits::GridViewType::dimension };
820 
821  typedef typename T::Traits::RangeFieldType Real;
822  typedef typename ConvectionDiffusionBoundaryConditions::Type BCType;
823 
824  public:
825  // pattern assembly flags
826  enum { doPatternVolume = false };
827  enum { doPatternSkeleton = false };
828 
829  // residual assembly flags
830  enum { doAlphaVolume = true };
831  enum { doAlphaBoundary = true };
832 
834  // supply time step from implicit Euler scheme
835  ConvectionDiffusionTemporalResidualEstimator (T& param_, double time_, double dt_)
836  : param(param_), time(time_), dt(dt_)
837  {}
838 
839  // volume integral depending on test and ansatz functions
840  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
841  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
842  {
843  // domain and range field type
844  typedef typename LFSU::Traits::FiniteElementType::
845  Traits::LocalBasisType::Traits::DomainFieldType DF;
846  typedef typename LFSU::Traits::FiniteElementType::
847  Traits::LocalBasisType::Traits::RangeFieldType RF;
848  typedef typename LFSU::Traits::FiniteElementType::
849  Traits::LocalBasisType::Traits::RangeType RangeType;
850  typedef typename LFSU::Traits::SizeType size_type;
851  typedef typename LFSU::Traits::FiniteElementType::
852  Traits::LocalBasisType::Traits::JacobianType JacobianType;
853 
854  // dimensions
855  const int dim = EG::Geometry::dimension;
856  const int intorder = 2*lfsu.finiteElement().localBasis().order();
857 
858  // select quadrature rule
859  Dune::GeometryType gt = eg.geometry().type();
860  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
861 
862  // interpolate f as finite element function to compute the gradient
863  CD_RHS_LocalAdapter<T,EG> f_adapter(param,eg);
864  std::vector<RF> f_up, f_down, f_mid;
865  param.setTime(time);
866  lfsu.finiteElement().localInterpolation().interpolate(f_adapter,f_down);
867  param.setTime(time+0.5*dt);
868  lfsu.finiteElement().localInterpolation().interpolate(f_adapter,f_mid);
869  param.setTime(time+dt);
870  lfsu.finiteElement().localInterpolation().interpolate(f_adapter,f_up);
871 
872  // loop over quadrature points
873  RF sum(0.0);
874  RF sum_grad(0.0);
875  RF fsum_grad_up(0.0);
876  RF fsum_grad_mid(0.0);
877  RF fsum_grad_down(0.0);
878  for (const auto& ip : rule)
879  {
880  // evaluate basis functions
881  std::vector<RangeType> phi(lfsu.size());
882  lfsu.finiteElement().localBasis().evaluateFunction(ip.position(),phi);
883 
884  // evaluate u
885  RF u=0.0;
886  for (size_type i=0; i<lfsu.size(); i++)
887  u += x(lfsu,i)*phi[i];
888 
889  // integrate jump
890  RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
891  sum += u*u*factor;
892 
893  // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
894  std::vector<JacobianType> js(lfsu.size());
895  lfsu.finiteElement().localBasis().evaluateJacobian(ip.position(),js);
896 
897  // transform gradients of shape functions to real element
898  const typename EG::Geometry::JacobianInverseTransposed jac =
899  eg.geometry().jacobianInverseTransposed(ip.position());
900  std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
901  for (size_type i=0; i<lfsu.size(); i++)
902  jac.mv(js[i][0],gradphi[i]);
903 
904  // compute gradient of u
905  Dune::FieldVector<RF,dim> gradu(0.0);
906  for (size_type i=0; i<lfsu.size(); i++)
907  gradu.axpy(x(lfsu,i),gradphi[i]);
908 
909  // integrate jump of gradient
910  sum_grad += (gradu*gradu)*factor;
911 
912  // compute gradients of f
913  Dune::FieldVector<RF,dim> gradf_down(0.0);
914  for (size_type i=0; i<lfsu.size(); i++) gradf_down.axpy(f_down[i],gradphi[i]);
915  Dune::FieldVector<RF,dim> gradf_mid(0.0);
916  for (size_type i=0; i<lfsu.size(); i++) gradf_mid.axpy(f_mid[i],gradphi[i]);
917  Dune::FieldVector<RF,dim> gradf_up(0.0);
918  for (size_type i=0; i<lfsu.size(); i++) gradf_up.axpy(f_up[i],gradphi[i]);
919  Dune::FieldVector<RF,dim> gradf_average(0.0);
920  for (size_type i=0; i<lfsu.size(); i++)
921  gradf_average.axpy((1.0/6.0)*f_down[i]+(2.0/3.0)*f_mid[i]+(1.0/6.0)*f_up[i],gradphi[i]);
922 
923  // integrate grad(f-f_average)
924  gradf_down -= gradf_average;
925  fsum_grad_down += (gradf_down*gradf_down)*factor;
926  gradf_mid -= gradf_average;
927  fsum_grad_mid += (gradf_mid*gradf_mid)*factor;
928  gradf_up -= gradf_average;
929  fsum_grad_up += (gradf_up*gradf_up)*factor;
930  }
931 
932  // accumulate cell indicator
933  DF h_T = diameter(eg.geometry());
934  r.accumulate(lfsv,0,dt * sum_grad); // k_n*||grad(jump)||^2
935  r.accumulate(lfsv,0,dt*dt * dt*((1.0/6.0)*fsum_grad_down+(2.0/3.0)*fsum_grad_mid+(1.0/6.0)*fsum_grad_up)); // k_n^2*||grad(f-time_average(f))||^2_s_t
936  }
937 
938  // boundary integral depending on test and ansatz functions
939  // We put the Dirchlet evaluation also in the alpha term to save some geometry evaluations
940  template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
941  void alpha_boundary (const IG& ig,
942  const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
943  R& r_s) const
944  {
945  // domain and range field type
946  typedef typename LFSU::Traits::FiniteElementType::
947  Traits::LocalBasisType::Traits::DomainFieldType DF;
948  typedef typename LFSU::Traits::FiniteElementType::
949  Traits::LocalBasisType::Traits::RangeFieldType RF;
950 
951  // dimensions
952  const int dim = IG::dimension;
953 
954  // select quadrature rule
955  const int intorder = 2*lfsu_s.finiteElement().localBasis().order();
956  Dune::GeometryType gtface = ig.geometryInInside().type();
957  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
958 
959  // evaluate boundary condition
960  const Dune::FieldVector<DF,dim-1>
961  face_local = Dune::ReferenceElements<DF,dim-1>::general(gtface).position(0,0);
962  BCType bctype = param.bctype(ig.intersection(),face_local);
964  return;
965 
966  // loop over quadrature points and integrate
967  RF sum_up(0.0);
968  RF sum_down(0.0);
969  for (const auto& ip : rule)
970  {
971  // evaluate flux boundary condition
972  param.setTime(time);
973  RF j_down = param.j(ig.intersection(),ip.position());
974  param.setTime(time+0.5*dt);
975  RF j_mid = param.j(ig.intersection(),ip.position());
976  param.setTime(time+dt);
977  RF j_up = param.j(ig.intersection(),ip.position());
978 
979  // integrate
980  RF factor = ip.weight() * ig.geometry().integrationElement(ip.position());
981  sum_down += (j_down-j_mid)*(j_down-j_mid)*factor;
982  sum_up += (j_up-j_mid)*(j_up-j_mid)*factor;
983  }
984 
985  // accumulate indicator
986  //DF h_T = diameter(ig.geometry());
987  auto inside_cell = ig.inside();
988  DF h_T = diameter(inside_cell.geometry());
989  r_s.accumulate(lfsv_s,0,(h_T+dt*dt/h_T)*dt*0.5*(sum_down+sum_up));
990  }
991 
992  private:
993  T& param; // two phase parameter class
994  double time;
995  double dt;
996 
997  template<class GEO>
998  typename GEO::ctype diameter (const GEO& geo) const
999  {
1000  typedef typename GEO::ctype DF;
1001  DF hmax = -1.0E00;
1002  const int dim = GEO::coorddimension;
1003  for (int i=0; i<geo.corners(); i++)
1004  {
1005  Dune::FieldVector<DF,dim> xi = geo.corner(i);
1006  for (int j=i+1; j<geo.corners(); j++)
1007  {
1008  Dune::FieldVector<DF,dim> xj = geo.corner(j);
1009  xj -= xi;
1010  hmax = std::max(hmax,xj.two_norm());
1011  }
1012  }
1013  return hmax;
1014  }
1015 
1016  };
1017 
1018 
1019 
1020  }
1021 }
1022 #endif
ConvectionDiffusionFEMResidualEstimator(T &param_)
constructor: pass parameter object
Definition: convectiondiffusionfem.hh:390
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: convectiondiffusionfem.hh:941
Definition: convectiondiffusionfem.hh:50
void setTime(double t)
set time in parameter class
Definition: convectiondiffusionfem.hh:341
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: convectiondiffusionfem.hh:540
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: convectiondiffusionfem.hh:203
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
Definition: convectiondiffusionfem.hh:53
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusionfem.hh:138
const IG & ig
Definition: constraints.hh:147
Type
Definition: convectiondiffusionparameter.hh:67
ConvectionDiffusionTemporalResidualEstimator1(T &param_, double time_, double dt_)
constructor: pass parameter object
Definition: convectiondiffusionfem.hh:677
Default flags for all local operators.
Definition: flags.hh:18
void clearCmax()
Definition: convectiondiffusionfem.hh:743
Definition: convectiondiffusionparameter.hh:67
Definition: convectiondiffusionfem.hh:658
sparsity pattern generator
Definition: pattern.hh:13
Definition: convectiondiffusionparameter.hh:67
ConvectionDiffusionTemporalResidualEstimator(T &param_, double time_, double dt_)
constructor: pass parameter object
Definition: convectiondiffusionfem.hh:835
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionfem.hh:63
ConvectionDiffusionFEM(T &param_, int intorderadd_=0)
Definition: convectiondiffusionfem.hh:56
Definition: convectiondiffusionparameter.hh:67
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: convectiondiffusionfem.hh:448
Definition: convectiondiffusionfem.hh:39
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionfem.hh:396
void jacobian_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, M &mat_s) const
Definition: convectiondiffusionfem.hh:282
Definition: adaptivity.hh:27
Definition: convectiondiffusionfem.hh:815
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
Implement jacobian_boundary() based on alpha_boundary()
Definition: defaultimp.hh:251
double getCmax() const
Definition: convectiondiffusionfem.hh:748
static const int dim
Definition: adaptivity.hh:83
Definition: convectiondiffusionfem.hh:371
Definition: convectiondiffusionfem.hh:782
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:15
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionfem.hh:841
Definition: convectiondiffusionfem.hh:54
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusionfem.hh:683
void evaluate(const X &x, Y &y) const
Definition: convectiondiffusionfem.hh:789
CD_RHS_LocalAdapter(const T &t_, const EG &eg_)
Definition: convectiondiffusionfem.hh:785
const EG & eg
Definition: constraints.hh:280
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538