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