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