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