dune-pdelab 2.7-git
Loading...
Searching...
No Matches
nonlinearconvectiondiffusionfem.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil -*-
2#ifndef DUNE_PDELAB_LOCALOPERATOR_NONLINEARCONVECTIONDIFFUSIONFEM_HH
3#define DUNE_PDELAB_LOCALOPERATOR_NONLINEARCONVECTIONDIFFUSIONFEM_HH
4
5#include<vector>
6
7#include<dune/common/exceptions.hh>
8#include<dune/common/fvector.hh>
9#include<dune/geometry/type.hh>
10
11#include<dune/geometry/referenceelements.hh>
12
16
17#include"defaultimp.hh"
18#include"pattern.hh"
19#include"flags.hh"
20#include"idefault.hh"
21
22
23namespace Dune {
24 namespace PDELab {
28
39 template<typename GV, typename RF>
41 {
43 using GridViewType = GV;
44
46 enum {
48 dimDomain = GV::dimension
49 };
50
52 using DomainFieldType = typename GV::Grid::ctype;
53
55 using DomainType = Dune::FieldVector<DomainFieldType,dimDomain>;
56
58 using IntersectionDomainType = Dune::FieldVector<DomainFieldType,dimDomain-1>;
59
61 using RangeFieldType = RF;
62
64 using RangeType = Dune::FieldVector<RF,GV::dimensionworld>;
65
67 using PermTensorType = Dune::FieldMatrix<RangeFieldType,dimDomain,dimDomain>;
68
70 using ElementType = typename GV::Traits::template Codim<0>::Entity;
71 using IntersectionType = typename GV::Intersection;
72 };
73
75 template<class T, class Imp>
77 {
78 public:
79 using Traits = T;
80
82 typename Traits::RangeFieldType
83 f (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
84 typename Traits::RangeFieldType u) const
85 {
86 return asImp().f(e,x,u);
87 }
88
90 typename Traits::RangeFieldType
91 w (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
92 typename Traits::RangeFieldType u) const
93 {
94 return asImp().w(e,x,u);
95 }
96
98 typename Traits::RangeFieldType
99 v (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
100 typename Traits::RangeFieldType u) const
101 {
102 return asImp().v(e,x,u);
103 }
104
106 typename Traits::PermTensorType
107 D (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
108 {
109 return asImp().D(e,x);
110 }
111
113 typename Traits::RangeType
114 q (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
115 typename Traits::RangeFieldType u) const
116 {
117 return asImp().q(e,x,u);
118 }
119
120 template<typename I>
122 const I & intersection, /*@\label{bcp:name}@*/
123 const Dune::FieldVector<typename I::ctype, I::mydimension> & coord
124 ) const
125 {
126 return asImp().isDirichlet( intersection, coord );
127 }
128
130 typename Traits::RangeFieldType
131 g (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
132 {
133 return asImp().g(e,x);
134 }
135
137 // Good: The dependence on u allows us to implement Robin type boundary conditions.
138 // Bad: This interface cannot be used for mixed finite elements where the flux is the essential b.c.
139 typename Traits::RangeFieldType
140 j (const typename Traits::ElementType& e, const typename Traits::DomainType& x,
141 typename Traits::RangeFieldType u) const
142 {
143 return asImp().j(e,x);
144 }
145
146 private:
147 Imp& asImp () {return static_cast<Imp &> (*this);}
148 const Imp& asImp () const {return static_cast<const Imp &>(*this);}
149 };
150
151
156 template<typename T>
158 : public Dune::PDELab::DirichletConstraintsParameters /*@\label{bcp:base}@*/
159 {
160 const typename T::Traits::GridViewType gv;
161 const T& t;
162
163 public:
164 BCTypeParam_CD( const typename T::Traits::GridViewType& gv_, const T& t_ )
165 : gv( gv_ ), t( t_ )
166 {
167 }
168
169 template<typename I>
171 const I & intersection, /*@\label{bcp:name}@*/
172 const Dune::FieldVector<typename I::ctype, I::mydimension> & coord
173 ) const
174 {
175 return t.isDirichlet( intersection, coord );
176 }
177 };
178
179
184 template<typename T>
186 : public GridFunctionBase<Dune::PDELab::GridFunctionTraits<typename T::Traits::GridViewType,
187 typename T::Traits::RangeFieldType,
188 1,Dune::FieldVector<typename T::Traits::RangeFieldType,1> >
189 ,DirichletBoundaryCondition_CD<T> >
190 {
191 public:
192 using Traits = Dune::PDELab::GridFunctionTraits<typename T::Traits::GridViewType,
193 typename T::Traits::RangeFieldType,
194 1,Dune::FieldVector<typename T::Traits::RangeFieldType,1> >;
195
197 DirichletBoundaryCondition_CD (const typename Traits::GridViewType& g_, const T& t_) : g(g_), t(t_) {}
198
200 inline void evaluate (const typename Traits::ElementType& e,
201 const typename Traits::DomainType& x,
202 typename Traits::RangeType& y) const
203 {
204 y = t.g(e,x);
205 }
206
207 inline const typename Traits::GridViewType& getGridView () const
208 {
209 return g;
210 }
211
212 private:
213 typename Traits::GridViewType g;
214 const T& t;
215 };
216
217
223 template<typename T>
225 public NumericalJacobianApplyVolume<NonLinearConvectionDiffusionFEM<T> >,
226 public NumericalJacobianApplyBoundary<NonLinearConvectionDiffusionFEM<T> >,
227 public NumericalJacobianVolume<NonLinearConvectionDiffusionFEM<T> >,
228 public NumericalJacobianBoundary<NonLinearConvectionDiffusionFEM<T> >,
229 public FullVolumePattern,
232 {
233 public:
234 // pattern assembly flags
235 enum { doPatternVolume = true };
236
237 // residual assembly flags
238 enum { doAlphaVolume = true };
239 enum { doAlphaBoundary = true };
240
241 NonLinearConvectionDiffusionFEM (T& param_, int intorder_=2)
242 : param(param_), intorder(intorder_)
243 {}
244
245 // volume integral depending on test and ansatz functions
246 template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
247 void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
248 {
249 // define types
250 using RF = typename LFSU::Traits::FiniteElementType::
251 Traits::LocalBasisType::Traits::RangeFieldType;
252 using JacobianType = typename LFSU::Traits::FiniteElementType::
253 Traits::LocalBasisType::Traits::JacobianType;
254 using RangeType = typename LFSU::Traits::FiniteElementType::
255 Traits::LocalBasisType::Traits::RangeType;
256 using size_type = typename LFSU::Traits::SizeType;
257
258 // dimensions
259 const int dim = EG::Geometry::mydimension;
260
261 // Reference to cell
262 const auto& cell = eg.entity();
263
264 // select quadrature rule
265 auto geo = eg.geometry();
266
267 // evaluate diffusion tensor at cell center, assume it is constant over elements
268 auto ref_el = referenceElement(geo);
269 auto localcenter = ref_el.position(0,0);
270 auto tensor = param.D(cell,localcenter);
271
272 // evaluate nonlinearity w(x_i); we assume here it is a Lagrange basis!
273 std::vector<typename T::Traits::RangeFieldType> w(lfsu.size());
274 for (size_type i=0; i<lfsu.size(); i++)
275 w[i] = param.w(cell,localcenter,x(lfsu,i));
276
277 // Transformation
278 typename EG::Geometry::JacobianInverseTransposed jac;
279
280 // Initialize vectors outside for loop
281 std::vector<RangeType> phi(lfsu.size());
282 std::vector<JacobianType> js(lfsu.size());
283 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
284 Dune::FieldVector<RF,dim> vgradu(0.0);
285 Dune::FieldVector<RF,dim> Dvgradu(0.0);
286
287 // loop over quadrature points
288 for (const auto& ip : quadratureRule(geo,intorder))
289 {
290 // evaluate basis functions
291 lfsu.finiteElement().localBasis().evaluateFunction(ip.position(),phi);
292
293 // evaluate u
294 RF u=0.0;
295 for (size_type i=0; i<lfsu.size(); i++)
296 u += w[i]*phi[i];
297
298 // evaluate source term
299 auto f = param.f(cell,ip.position(),u);
300
301 // evaluate flux term
302 auto q = param.q(cell,ip.position(),u);
303
304 // evaluate gradient of shape functions (we assume Galerkin method lfsu=lfsv)
305 lfsu.finiteElement().localBasis().evaluateJacobian(ip.position(),js);
306
307 // transform gradients of shape functions to real element
308 jac = geo.jacobianInverseTransposed(ip.position());
309 for (size_type i=0; i<lfsu.size(); i++)
310 {
311 gradphi[i] = 0.0;
312 jac.umv(js[i][0],gradphi[i]);
313 }
314
315 // v(u) compute gradient of u
316 vgradu = 0.0;
317 for (size_type i=0; i<lfsu.size(); i++)
318 vgradu.axpy(w[i],gradphi[i]);
319 vgradu *= param.v(cell,ip.position(),u);
320
321 // compute D * v(u) * gradient of u
322 Dvgradu = 0.0;
323 tensor.umv(vgradu,Dvgradu);
324
325 // integrate (K grad u)*grad phi_i + a_0*u*phi_i
326 auto factor = ip.weight() * geo.integrationElement(ip.position());
327 for (size_type i=0; i<lfsu.size(); i++)
328 r.accumulate(lfsu,i,( Dvgradu*gradphi[i] - q*gradphi[i] - f*phi[i] )*factor);
329 }
330 }
331
332 // boundary integral
333 template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
334 void alpha_boundary (const IG& ig,
335 const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
336 R& r_s) const
337 {
338 // define types
339 using RF = typename LFSV::Traits::FiniteElementType::
340 Traits::LocalBasisType::Traits::RangeFieldType;
341 using RangeType = typename LFSV::Traits::FiniteElementType::
342 Traits::LocalBasisType::Traits::RangeType;
343 using size_type = typename LFSV::Traits::SizeType;
344
345 // get inside cell entity
346 const auto& cell_inside = ig.inside();
347
348 // get geometry
349 auto geo = ig.geometry();
350
351 // Get geometry of intersection in local coordinates of cell_inside
352 auto geo_in_inside = ig.geometryInInside();
353
354 // evaluate nonlinearity w(x_i); we assume here it is a Lagrange basis!
355 auto ref_el_in_inside = referenceElement(geo_in_inside);
356 auto local_face_center = ref_el_in_inside.position(0,0);
357 auto face_center_in_element = geo_in_inside.global(local_face_center);
358 std::vector<typename T::Traits::RangeFieldType> w(lfsu_s.size());
359 for (size_type i=0; i<lfsu_s.size(); i++)
360 w[i] = param.w(cell_inside,face_center_in_element,x_s(lfsu_s,i));
361
362 // Initialize vectors outside for loop
363 std::vector<RangeType> phi(lfsv_s.size());
364
365 // loop over quadrature points and integrate normal flux
366 for (const auto& ip : quadratureRule(geo,intorder))
367 {
368 // evaluate boundary condition type
369 // skip rest if we are on Dirichlet boundary
370 if( param.isDirichlet( ig.intersection(), ip.position() ) )
371 continue;
372
373 // position of quadrature point in local coordinates of element
374 auto local = geo_in_inside.global(ip.position());
375
376 // evaluate test shape functions
377 lfsv_s.finiteElement().localBasis().evaluateFunction(local,phi);
378
379 // evaluate u
380 RF u=0.0;
381 for (size_type i=0; i<lfsu_s.size(); i++)
382 u += w[i]*phi[i];
383
384 // evaluate flux boundary condition
385 auto j = param.j(cell_inside,local,u);
386
387 // integrate j
388 auto factor = ip.weight()*geo.integrationElement(ip.position());
389 for (size_type i=0; i<lfsv_s.size(); i++)
390 r_s.accumulate(lfsu_s,i,j*phi[i]*factor);
391 }
392 }
393
395 void setTime (double t)
396 {
397 param.setTime(t);
398 }
399
400 private:
401 T& param;
402 int intorder;
403 };
404
406 } // namespace PDELab
407} // namespace Dune
408
409#endif // DUNE_PDELAB_LOCALOPERATOR_NONLINEARCONVECTIONDIFFUSIONFEM_HH
VTKWriter & w
Definition: function.hh:842
static const int dim
Definition: adaptivity.hh:84
const IG & ig
Definition: constraints.hh:149
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
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::FieldVector< GV::Grid::ctype, GV::dimension > DomainType
domain type in dim-size coordinates
Definition: function.hh:50
GV::Traits::template Codim< 0 >::Entity ElementType
codim 0 entity
Definition: function.hh:119
GV GridViewType
The type of the grid view the function lives on.
Definition: function.hh:116
traits class holding the function signature, same as in local function
Definition: function.hh:183
leaf of a function tree
Definition: function.hh:302
Definition: constraintsparameters.hh:26
Default flags for all local operators.
Definition: flags.hh:19
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
traits class for two phase parameter class
Definition: nonlinearconvectiondiffusionfem.hh:41
Dune::FieldVector< DomainFieldType, dimDomain > DomainType
domain type
Definition: nonlinearconvectiondiffusionfem.hh:55
typename GV::Grid::ctype DomainFieldType
Export type for domain field.
Definition: nonlinearconvectiondiffusionfem.hh:52
Dune::FieldVector< DomainFieldType, dimDomain-1 > IntersectionDomainType
domain type
Definition: nonlinearconvectiondiffusionfem.hh:58
GV GridViewType
the grid view
Definition: nonlinearconvectiondiffusionfem.hh:43
@ dimDomain
dimension of the domain
Definition: nonlinearconvectiondiffusionfem.hh:48
typename GV::Traits::template Codim< 0 >::Entity ElementType
grid types
Definition: nonlinearconvectiondiffusionfem.hh:70
RF RangeFieldType
Export type for range field.
Definition: nonlinearconvectiondiffusionfem.hh:61
Dune::FieldVector< RF, GV::dimensionworld > RangeType
range type
Definition: nonlinearconvectiondiffusionfem.hh:64
Dune::FieldMatrix< RangeFieldType, dimDomain, dimDomain > PermTensorType
permeability tensor type
Definition: nonlinearconvectiondiffusionfem.hh:67
typename GV::Intersection IntersectionType
Definition: nonlinearconvectiondiffusionfem.hh:71
base class for parameter class
Definition: nonlinearconvectiondiffusionfem.hh:77
bool isDirichlet(const I &intersection, const Dune::FieldVector< typename I::ctype, I::mydimension > &coord) const
Definition: nonlinearconvectiondiffusionfem.hh:121
Traits::RangeFieldType g(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
Dirichlet boundary condition value.
Definition: nonlinearconvectiondiffusionfem.hh:131
Traits::RangeFieldType j(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType u) const
Neumann boundary condition.
Definition: nonlinearconvectiondiffusionfem.hh:140
Traits::RangeFieldType v(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType u) const
scalar nonlinearity in diffusion coefficient
Definition: nonlinearconvectiondiffusionfem.hh:99
Traits::RangeFieldType f(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType u) const
source/reaction term
Definition: nonlinearconvectiondiffusionfem.hh:83
T Traits
Definition: nonlinearconvectiondiffusionfem.hh:79
Traits::RangeFieldType w(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType u) const
nonlinearity under gradient
Definition: nonlinearconvectiondiffusionfem.hh:91
Traits::RangeType q(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeFieldType u) const
nonlinear flux vector
Definition: nonlinearconvectiondiffusionfem.hh:114
Traits::PermTensorType D(const typename Traits::ElementType &e, const typename Traits::DomainType &x) const
tensor diffusion coefficient
Definition: nonlinearconvectiondiffusionfem.hh:107
Definition: nonlinearconvectiondiffusionfem.hh:159
BCTypeParam_CD(const typename T::Traits::GridViewType &gv_, const T &t_)
Definition: nonlinearconvectiondiffusionfem.hh:164
bool isDirichlet(const I &intersection, const Dune::FieldVector< typename I::ctype, I::mydimension > &coord) const
Definition: nonlinearconvectiondiffusionfem.hh:170
Definition: nonlinearconvectiondiffusionfem.hh:190
DirichletBoundaryCondition_CD(const typename Traits::GridViewType &g_, const T &t_)
constructor
Definition: nonlinearconvectiondiffusionfem.hh:197
const Traits::GridViewType & getGridView() const
Definition: nonlinearconvectiondiffusionfem.hh:207
void evaluate(const typename Traits::ElementType &e, const typename Traits::DomainType &x, typename Traits::RangeType &y) const
Evaluate the GridFunction at given position.
Definition: nonlinearconvectiondiffusionfem.hh:200
Definition: nonlinearconvectiondiffusionfem.hh:232
NonLinearConvectionDiffusionFEM(T &param_, int intorder_=2)
Definition: nonlinearconvectiondiffusionfem.hh:241
@ doPatternVolume
Definition: nonlinearconvectiondiffusionfem.hh:235
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: nonlinearconvectiondiffusionfem.hh:334
void setTime(double t)
set time in parameter class
Definition: nonlinearconvectiondiffusionfem.hh:395
@ doAlphaBoundary
Definition: nonlinearconvectiondiffusionfem.hh:239
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: nonlinearconvectiondiffusionfem.hh:247
@ doAlphaVolume
Definition: nonlinearconvectiondiffusionfem.hh:238
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:32
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:251
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume()
Definition: numericaljacobianapply.hh:34
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary()
Definition: numericaljacobianapply.hh:287
sparsity pattern generator
Definition: pattern.hh:14