DUNE PDELab (git)

diffusionmixed.hh
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
11
13#include<dune/geometry/referenceelements.hh>
14
15#include<dune/pdelab/common/quadraturerules.hh>
16
17#include"defaultimp.hh"
18#include"pattern.hh"
19#include"flags.hh"
20#include "convectiondiffusionparameter.hh"
21
22namespace 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,
40 public LocalOperatorDefaultFlags
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
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
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
275 if (bctype == ConvectionDiffusionBoundaryConditions::Neumann)
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
vector space out of a tensor product of fields.
Definition: fvector.hh:91
A few common exception classes.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
constexpr index_constant< 0 > _0
Compile time index with value 0.
Definition: indices.hh:52
constexpr index_constant< 1 > _1
Compile time index with value 1.
Definition: indices.hh:55
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
ImplementationDefined child(Node &&node, Indices... indices)
Extracts the child of a node given by a sequence of compile-time and run-time indices.
Definition: childextraction.hh:128
Namespace with predefined compile time indices for the range [0,19].
Definition: indices.hh:50
Dune namespace.
Definition: alignedallocator.hh:13
A unique label for each type of element that can occur in a grid.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)