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