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