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