DUNE PDELab (git)

nonlinearconvectiondiffusionfem.hh
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
10
11#include<dune/geometry/referenceelements.hh>
12
13#include<dune/pdelab/common/quadraturerules.hh>
14#include <dune/pdelab/common/function.hh>
15#include<dune/pdelab/constraints/common/constraintsparameters.hh>
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
56
59
61 using RangeFieldType = RF;
62
65
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>
121 bool isDirichlet(
122 const I & intersection, /*@\label{bcp:name}@*/
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>
170 bool isDirichlet(
171 const I & intersection, /*@\label{bcp:name}@*/
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,
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
derived_type & axpy(const field_type &a, const DenseVector< Other > &x)
vector space axpy operation ( *this += a x )
Definition: densevector.hh:575
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Definition: nonlinearconvectiondiffusionfem.hh:159
Definition: nonlinearconvectiondiffusionfem.hh:190
DirichletBoundaryCondition_CD(const typename Traits::GridViewType &g_, const T &t_)
constructor
Definition: nonlinearconvectiondiffusionfem.hh:197
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
sparsity pattern generator
Definition: pattern.hh:14
leaf of a function tree
Definition: function.hh:302
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
Default flags for all local operators.
Definition: flags.hh:19
Definition: nonlinearconvectiondiffusionfem.hh:232
void setTime(double t)
set time in parameter class
Definition: nonlinearconvectiondiffusionfem.hh:395
base class for parameter class
Definition: nonlinearconvectiondiffusionfem.hh:77
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
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
Implements linear and nonlinear versions of jacobian_apply_boundary() based on alpha_boundary()
Definition: numericaljacobianapply.hh:287
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume()
Definition: numericaljacobianapply.hh:34
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:251
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:32
A few common exception classes.
Implements a vector constructed from a given type representing a field and a compile-time given size.
unspecified value type referenceElement(T &&... t)
Returns a reference element for the objects t....
Dune namespace.
Definition: alignedallocator.hh:13
Static tag representing a codimension.
Definition: dimension.hh:24
Definition: constraintsparameters.hh:26
traits class holding the function signature, same as in local function
Definition: function.hh:183
traits class for two phase parameter class
Definition: nonlinearconvectiondiffusionfem.hh:41
typename GV::Grid::ctype DomainFieldType
Export type for domain field.
Definition: nonlinearconvectiondiffusionfem.hh:52
@ dimDomain
dimension of the domain
Definition: nonlinearconvectiondiffusionfem.hh:48
GV GridViewType
the grid view
Definition: nonlinearconvectiondiffusionfem.hh:43
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
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
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)