dune-pdelab  2.4-dev
diffusionmixed.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_PDELAB_DIFFUSIONMIXED_HH
3 #define DUNE_PDELAB_DIFFUSIONMIXED_HH
4 
5 #include <cstddef>
6 #include<vector>
7 
8 #include<dune/common/exceptions.hh>
9 #include<dune/common/fvector.hh>
10 #include<dune/common/fmatrix.hh>
11 
12 #include<dune/geometry/type.hh>
13 #include<dune/geometry/quadraturerules.hh>
14 #include<dune/geometry/referenceelements.hh>
15 
16 #include"defaultimp.hh"
17 #include"pattern.hh"
18 #include"flags.hh"
20 
21 namespace Dune {
22  namespace PDELab {
23 
24  // a local operator for solving the Poisson equation
25  // div sigma +a_0 u = f in \Omega,
26  // sigma = -K grad u in \Omega,
27  // u = g on \partial\Omega_D
28  // sigma \cdot \nu = j on \partial\Omega_N
29  // with H(div) conforming (mixed) finite elements
30  // param.A : diffusion tensor dependent on position
31  // param.c : Helmholtz term
32  // param.f : grid function type giving f
33  // param.bctype : grid function type selecting boundary condition
34  // param.g : grid function type giving g
35  template<typename PARAM>
36  class DiffusionMixed : public NumericalJacobianApplyVolume<DiffusionMixed<PARAM> >,
37  public NumericalJacobianVolume<DiffusionMixed<PARAM> >,
38  public FullVolumePattern,
40  {
41 
43 
44  public:
45  // pattern assembly flags
46  enum { doPatternVolume = true };
47 
48  // residual assembly flags
49  enum { doAlphaVolume = true };
50  enum { doLambdaVolume = true };
51  enum { doLambdaBoundary = true };
52 
53  DiffusionMixed ( const PARAM& param_,
54  int qorder_v_=2,
55  int qorder_p_=1 )
56  : param(param_),
57  qorder_v(qorder_v_),
58  qorder_p(qorder_p_)
59  {
60  }
61 
62  // volume integral depending on test and ansatz functions
63  template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
64  void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
65  {
66  // select the two components
67  typedef typename LFSU::template Child<0>::Type VelocitySpace;
68  const VelocitySpace& velocityspace = lfsu.template child<0>();
69  typedef typename LFSU::template Child<1>::Type PressureSpace;
70  const PressureSpace& pressurespace = lfsu.template child<1>();
71 
72  // domain and range field type
73  typedef typename VelocitySpace::Traits::FiniteElementType::
74  Traits::LocalBasisType::Traits::DomainFieldType DF;
75  typedef typename VelocitySpace::Traits::FiniteElementType::
76  Traits::LocalBasisType::Traits::RangeFieldType RF;
77  typedef typename VelocitySpace::Traits::FiniteElementType::
78  Traits::LocalBasisType::Traits::JacobianType VelocityJacobianType;
79  typedef typename VelocitySpace::Traits::FiniteElementType::
80  Traits::LocalBasisType::Traits::RangeType VelocityRangeType;
81  typedef typename PressureSpace::Traits::FiniteElementType::
82  Traits::LocalBasisType::Traits::RangeType PressureRangeType;
83 
84  // dimensions
85  const int dim = EG::Geometry::dimension;
86 
87  // evaluate transformation which must be linear
88  Dune::FieldVector<DF,dim> pos; pos = 0.0;
89  typename EG::Geometry::JacobianInverseTransposed jac;
90  jac = eg.geometry().jacobianInverseTransposed(pos);
91  jac.invert();
92  RF det = eg.geometry().integrationElement(pos);
93 
94  // evaluate diffusion tensor at cell center, assume it is constant over elements
95  //typename K::Traits::RangeType tensor;
96  Dune::GeometryType gt = eg.geometry().type();
97  Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
98 
99  typename PARAM::Traits::PermTensorType tensor;
100  tensor = param.A(eg.entity(),localcenter);
101  tensor.invert(); // need iverse for mixed method
102 
103  // \sigma\cdot v term
104  const Dune::QuadratureRule<DF,dim>& vrule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder_v);
105  // loop over quadrature points
106  for (const auto& ip : vrule)
107  {
108  // evaluate shape functions at ip (this is a Galerkin method)
109  std::vector<VelocityRangeType> vbasis(velocityspace.size());
110  velocityspace.finiteElement().localBasis().evaluateFunction(ip.position(),vbasis);
111 
112  // transform basis vectors
113  std::vector<VelocityRangeType> vtransformedbasis(velocityspace.size());
114  for (std::size_t i=0; i<velocityspace.size(); i++)
115  {
116  vtransformedbasis[i] = 0.0;
117  jac.umtv(vbasis[i],vtransformedbasis[i]);
118  }
119 
120  // compute sigma
121  VelocityRangeType sigma; sigma = 0.0;
122  for (std::size_t i=0; i<velocityspace.size(); i++)
123  sigma.axpy(x(velocityspace,i),vtransformedbasis[i]);
124 
125  // K^{-1} * sigma
126  VelocityRangeType Kinvsigma;
127  tensor.mv(sigma,Kinvsigma);
128 
129  // integrate (K^{-1}*sigma) * phi_i
130  RF factor = ip.weight() / det;
131  for (std::size_t i=0; i<velocityspace.size(); i++)
132  r.accumulate(velocityspace,i,(Kinvsigma*vtransformedbasis[i])*factor);
133  }
134 
135  // u div v term, div sigma q term, a0*u term
136  const Dune::QuadratureRule<DF,dim>& prule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder_p);
137  // loop over quadrature points
138  for (const auto& ip : prule)
139  {
140  // evaluate shape functions at ip (this is a Galerkin method)
141  std::vector<VelocityJacobianType> vbasis(velocityspace.size());
142  velocityspace.finiteElement().localBasis().evaluateJacobian(ip.position(),vbasis);
143  std::vector<PressureRangeType> pbasis(pressurespace.size());
144  pressurespace.finiteElement().localBasis().evaluateFunction(ip.position(),pbasis);
145 
146  // compute u
147  PressureRangeType u; u = 0.0;
148  for (std::size_t i=0; i<pressurespace.size(); i++)
149  u.axpy(x(pressurespace,i),pbasis[i]);
150 
151  // evaluate Helmholtz term (reaction term)
152  typename PARAM::Traits::RangeFieldType a0value = param.c(eg.entity(),ip.position());
153 
154  // integrate a0 * u * q
155  RF factor = ip.weight();
156  for (std::size_t i=0; i<pressurespace.size(); i++)
157  r.accumulate(pressurespace,i,-a0value*u*pbasis[i]*factor);
158 
159  // compute divergence of velocity basis functions on reference element
160  std::vector<RF> divergence(velocityspace.size(),0.0);
161  for (std::size_t i=0; i<velocityspace.size(); i++)
162  for (int j=0; j<dim; j++)
163  divergence[i] += vbasis[i][j][j];
164 
165  // integrate sigma * phi_i
166  for (std::size_t i=0; i<velocityspace.size(); i++)
167  r.accumulate(velocityspace,i,-u*divergence[i]*factor);
168 
169  // compute divergence of sigma
170  RF divergencesigma = 0.0;
171  for (std::size_t i=0; i<velocityspace.size(); i++)
172  divergencesigma += x(velocityspace,i)*divergence[i];
173 
174  // integrate div sigma * q
175  for (std::size_t i=0; i<pressurespace.size(); i++)
176  r.accumulate(pressurespace,i,-divergencesigma*pbasis[i]*factor);
177  }
178  }
179 
180  // volume integral depending only on test functions
181  template<typename EG, typename LFSV, typename R>
182  void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
183  {
184  // select the two components
185  typedef typename LFSV::template Child<1>::Type PressureSpace;
186  const PressureSpace& pressurespace = lfsv.template child<1>();
187 
188  // domain and range field type
189  typedef typename PressureSpace::Traits::FiniteElementType::
190  Traits::LocalBasisType::Traits::DomainFieldType DF;
191  typedef typename PressureSpace::Traits::FiniteElementType::
192  Traits::LocalBasisType::Traits::RangeFieldType RF;
193  typedef typename PressureSpace::Traits::FiniteElementType::
194  Traits::LocalBasisType::Traits::RangeType PressureRangeType;
195 
196  // dimensions
197  const int dim = EG::Geometry::dimension;
198 
199  // select quadrature rule
200  Dune::GeometryType gt = eg.geometry().type();
201  const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder_p);
202  // loop over quadrature points
203  for (const auto& ip : rule)
204  {
205  // evaluate shape functions
206  std::vector<PressureRangeType> pbasis(pressurespace.size());
207  pressurespace.finiteElement().localBasis().evaluateFunction(ip.position(),pbasis);
208 
209  // evaluate right hand side parameter function
210  RF y = param.f(eg.entity(),ip.position());
211 
212  // integrate f
213  RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
214  for (std::size_t i=0; i<pressurespace.size(); i++)
215  r.accumulate(pressurespace,i,y*pbasis[i]*factor);
216  }
217  }
218 
219  // boundary integral independen of ansatz functions
220  template<typename IG, typename LFSV, typename R>
221  void lambda_boundary (const IG& ig, const LFSV& lfsv, R& r) const
222  {
223  // select the two components
224  typedef typename LFSV::template Child<0>::Type VelocitySpace;
225  const VelocitySpace& velocityspace = lfsv.template child<0>();
226 
227  // domain and range field type
228  typedef typename VelocitySpace::Traits::FiniteElementType::
229  Traits::LocalBasisType::Traits::DomainFieldType DF;
230  typedef typename VelocitySpace::Traits::FiniteElementType::
231  Traits::LocalBasisType::Traits::RangeFieldType RF;
232  typedef typename VelocitySpace::Traits::FiniteElementType::
233  Traits::LocalBasisType::Traits::RangeType VelocityRangeType;
234 
235  // dimensions
236  const int dim = IG::dimension;
237 
238  auto inside_cell = ig.inside();
239 
240  // evaluate transformation which must be linear
241  Dune::FieldVector<DF,dim> pos; pos = 0.0;
242  typename IG::Entity::Geometry::JacobianInverseTransposed jac;
243  jac = inside_cell.geometry().jacobianInverseTransposed(pos);
244  jac.invert();
245  RF det = inside_cell.geometry().integrationElement(pos);
246 
247  // select quadrature rule
248  Dune::GeometryType gtface = ig.geometryInInside().type();
249  const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder_v);
250 
251  // loop over quadrature points and integrate normal flux
252  for (const auto& ip : rule)
253  {
254  // evaluate boundary condition type
255  BCType bctype = param.bctype(ig.intersection(),ip.position());
256 
257  // skip rest if we are on Neumann boundary
259  continue;
260 
261  // position of quadrature point in local coordinates of element
262  Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
263 
264  // evaluate test shape functions
265  std::vector<VelocityRangeType> vbasis(velocityspace.size());
266  velocityspace.finiteElement().localBasis().evaluateFunction(local,vbasis);
267 
268  // transform basis vectors
269  std::vector<VelocityRangeType> vtransformedbasis(velocityspace.size());
270  for (std::size_t i=0; i<velocityspace.size(); i++)
271  {
272  vtransformedbasis[i] = 0.0;
273  jac.umtv(vbasis[i],vtransformedbasis[i]);
274  }
275 
276  // evaluate Dirichlet boundary condition
277  RF y = param.g(inside_cell,local);
278 
279  // integrate g v*normal
280  RF factor = ip.weight()*ig.geometry().integrationElement(ip.position())/det;
281  for (std::size_t i=0; i<velocityspace.size(); i++)
282  r.accumulate(velocityspace,i,y*(vtransformedbasis[i]*ig.unitOuterNormal(ip.position()))*factor);
283  }
284  }
285 
286  private:
287  const PARAM& param;
288  int qorder_v;
289  int qorder_p;
290  };
291 
293  } // namespace PDELab
294 } // namespace Dune
295 
296 #endif
Definition: diffusionmixed.hh:49
Definition: diffusionmixed.hh:50
Definition: diffusionmixed.hh:46
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
const IG & ig
Definition: constraints.hh:147
Type
Definition: convectiondiffusionparameter.hh:67
Default flags for all local operators.
Definition: flags.hh:18
sparsity pattern generator
Definition: pattern.hh:13
Definition: diffusionmixed.hh:51
Definition: convectiondiffusionparameter.hh:67
Definition: adaptivity.hh:27
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
static const int dim
Definition: adaptivity.hh:83
DiffusionMixed(const PARAM &param_, int qorder_v_=2, int qorder_p_=1)
Definition: diffusionmixed.hh:53
Definition: diffusionmixed.hh:36
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: diffusionmixed.hh:182
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusionmixed.hh:64
const EG & eg
Definition: constraints.hh:280
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: diffusionmixed.hh:221