DUNE PDELab (git)

dgnavierstokesvelvecfem.hh
1// -*- tab-width: 2; indent-tabs-mode: nil -*-
2// vi: set et ts=2 sw=2 sts=2:
3
4#ifndef DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKESVELVECFEM_HH
5#define DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKESVELVECFEM_HH
6
11
12#include <dune/localfunctions/common/interfaceswitch.hh>
13#include <dune/pdelab/localoperator/idefault.hh>
14
15#include <dune/pdelab/common/quadraturerules.hh>
16#include <dune/pdelab/localoperator/defaultimp.hh>
17#include <dune/pdelab/localoperator/pattern.hh>
18#include <dune/pdelab/localoperator/flags.hh>
19#include <dune/pdelab/localoperator/dgnavierstokesparameter.hh>
20#include <dune/pdelab/localoperator/navierstokesmass.hh>
21
22#ifndef VBLOCK
23#define VBLOCK 0
24#endif
25#define PBLOCK (- VBLOCK + 1)
26
27namespace Dune {
28 namespace PDELab {
29
30 template<class Basis, class Dummy = void>
31 struct VectorBasisInterfaceSwitch {
33 using DomainLocal = typename Basis::Traits::DomainLocal;
35 using RangeField = typename Basis::Traits::RangeField;
37 static const std::size_t dimRange = Basis::Traits::dimRange;
38
40 template<typename Geometry>
41 static void jacobian(const Basis& basis, const Geometry& geometry,
42 const DomainLocal& xl,
43 std::vector<FieldMatrix<RangeField, dimRange,
45 {
46 jac.resize(basis.size());
47 basis.evaluateJacobian(xl, jac);
48 }
49 };
50
52 template<class Basis>
53 struct VectorBasisInterfaceSwitch<
54 Basis, typename std::enable_if<
55 Dune::AlwaysTrue<
56 std::integral_constant<
57 std::size_t,
58 Basis::Traits::dimDomain
59 >
60 >::value
61 >::type
62 >
63 {
65 using DomainLocal = typename Basis::Traits::DomainType;
67 using RangeField = typename Basis::Traits::RangeFieldType;
69 static const std::size_t dimRange = Basis::Traits::dimRange;
70
72 template<typename Geometry>
73 static void jacobian(const Basis& basis, const Geometry& geometry,
74 const DomainLocal& xl,
75 std::vector<FieldMatrix<RangeField, dimRange,
77 {
78 jac.resize(basis.size());
79
80 std::vector<FieldMatrix<
81 RangeField, dimRange, Geometry::coorddimension> > ljac(basis.size());
82 basis.evaluateJacobian(xl, ljac);
83
84 const typename Geometry::JacobianInverseTransposed& geojac =
85 geometry.jacobianInverseTransposed(xl);
86
87 for(std::size_t i = 0; i < basis.size(); ++i)
88 for(std::size_t row=0; row < dimRange; ++row)
89 geojac.mv(ljac[i][row], jac[i][row]);
90 }
91 };
92
100 template<typename PRM>
104 public InstationaryLocalOperatorDefaultMethods<typename PRM::Traits::RangeField>
105 {
107 using RF = typename PRM::Traits::RangeField;
108
110 using Real = typename InstatBase::RealType;
111
112 static const bool navier = PRM::assemble_navier;
113 static const bool full_tensor = PRM::assemble_full_tensor;
114
115 public :
116
117 // pattern assembly flags
118 enum { doPatternVolume = true };
119 enum { doPatternSkeleton = true };
120
121 // call the assembler for each face only once
122 enum { doSkeletonTwoSided = false };
123
124 // residual assembly flags
125 enum { doAlphaVolume = true };
126 enum { doAlphaSkeleton = true };
127 enum { doAlphaBoundary = true };
128 enum { doLambdaVolume = true };
129
142 DGNavierStokesVelVecFEM (PRM& _prm, int _superintegration_order=0) :
143 prm(_prm), superintegration_order(_superintegration_order),
144 current_dt(1.0)
145 {}
146
147 // Store current dt
148 void preStep (Real , Real dt, int )
149 {
150 current_dt = dt;
151 }
152
153 // set time in parameter class
154 void setTime(Real t)
155 {
157 prm.setTime(t);
158 }
159
160 // volume integral depending on test and ansatz functions
161 template<typename EG, typename LFSU, typename X, typename LFSV, typename R>
162 void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const
163 {
164 const unsigned int dim = EG::Geometry::mydimension;
165
166 // subspaces
167 using namespace Indices;
168 using LFSV_V = TypeTree::Child<LFSV,_0>;
169 const auto& lfsv_v = child(lfsv,_0);
170 const auto& lfsu_v = child(lfsu,_0);
171
172 using LFSV_P = TypeTree::Child<LFSV,_1>;
173 const auto& lfsv_p = child(lfsv,_1);
174 const auto& lfsu_p = child(lfsu,_1);
175
176 // domain and range field type
177 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
178 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
179 using VectorBasisSwitch_V = VectorBasisInterfaceSwitch<typename FESwitch_V::Basis >;
180 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
181 using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
182 using RF = typename BasisSwitch_V::RangeField;
183 using Range_V = typename BasisSwitch_V::Range;
184 using Range_P = typename BasisSwitch_P::Range;
185 using size_type = typename LFSV::Traits::SizeType;
186
187 // Get geometry
188 auto geo = eg.geometry();
189
190 // Determine quadrature order
191 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
192 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
193 const int jac_order = geo.type().isSimplex() ? 0 : 1;
194 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
195
196 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
197
198 // loop over quadrature points
199 for (const auto& ip : quadratureRule(geo,qorder))
200 {
201 auto local = ip.position();
202 auto mu = prm.mu(eg,local);
203 auto rho = prm.rho(eg,local);
204
205 // compute u (if Navier term enabled)
206 std::vector<Range_V> phi_v(lfsv_v.size());
207 Range_V val_u(0.0);
208 if(navier) {
209 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
210 for (size_type i=0; i<lfsu_v.size(); i++)
211 val_u.axpy(x(lfsu_v,i),phi_v[i]);
212 }
213
214 // values of pressure shape functions
215 std::vector<Range_P> phi_p(lfsv_p.size());
216 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
217
218 // compute pressure value
219 Range_P val_p(0.0);
220 for (size_type i=0; i<lfsu_p.size(); i++)
221 val_p.axpy(x(lfsu_p,i),phi_p[i]);
222
223 // evaluate jacobian of velocity shape functions on reference element
224 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
225 VectorBasisSwitch_V::jacobian
226 (FESwitch_V::basis(lfsv_v.finiteElement()), geo, local, jac_phi_v);
227
228 // compute divergence of test functions
229 std::vector<RF> div_phi_v(lfsv_v.size(),0.0);
230 for (size_type i=0; i<lfsv_v.size(); i++)
231 for (size_type d=0; d<dim; d++)
232 div_phi_v[i] += jac_phi_v[i][d][d];
233
234 // compute velocity jacobian and divergence
236 RF div_u(0.0);
237 for (size_type i=0; i<lfsu_v.size(); i++){
238 jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
239 div_u += x(lfsu_v,i) * div_phi_v[i];
240 }
241
242 auto detj = geo.integrationElement(ip.position());
243 auto weight = ip.weight() * detj;
244
245 for (size_type i=0; i<lfsv_v.size(); i++) {
246 //================================================//
247 // \int (mu*grad_u*grad_v)
248 //================================================//
249 RF dvdu(0);
250 contraction(jac_u,jac_phi_v[i],dvdu);
251 r.accumulate(lfsv_v, i, dvdu * mu * weight);
252
253 //================================================//
254 // \int -p \nabla\cdot v
255 //================================================//
256 r.accumulate(lfsv_v, i, - div_phi_v[i] * val_p * weight);
257
258 //================================================//
259 // \int \rho ((u\cdot\nabla ) u )\cdot v
260 //================================================//
261 if(navier) {
262 // compute (grad u) u (matrix-vector product)
263 Range_V nabla_u_u(0.0);
264 jac_u.mv(val_u,nabla_u_u);
265 r.accumulate(lfsv_v, i, rho * (nabla_u_u*phi_v[i]) * weight);
266 } // end navier
267
268 } // end i
269
270 for (size_type i=0; i<lfsv_p.size(); i++) {
271 //================================================//
272 // \int -q \nabla\cdot u
273 //================================================//
274 r.accumulate(lfsv_p, i, - div_u * phi_p[i] * incomp_scaling * weight);
275 }
276
277 } // end loop quadrature points
278 } // end alpha_volume
279
280 // jacobian of volume term
281 template<typename EG, typename LFSU, typename X, typename LFSV,
282 typename LocalMatrix>
283 void jacobian_volume (const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv,
284 LocalMatrix& mat) const
285 {
286 const unsigned int dim = EG::Geometry::mydimension;
287
288 // subspaces
289 using namespace Indices;
290 using LFSV_V = TypeTree::Child<LFSV,_0>;
291 const auto& lfsv_v = child(lfsv,_0);
292 const auto& lfsu_v = child(lfsu,_0);
293
294 using LFSV_P = TypeTree::Child<LFSV,_1>;
295 const auto& lfsv_p = child(lfsv,_1);
296 const auto& lfsu_p = child(lfsu,_1);
297
298 // domain and range field type
299 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
300 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
301 using VectorBasisSwitch_V = VectorBasisInterfaceSwitch<typename FESwitch_V::Basis >;
302 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
303 using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
304 using RF = typename BasisSwitch_V::RangeField;
305 using Range_V = typename BasisSwitch_V::Range;
306 using Range_P = typename BasisSwitch_P::Range;
307 using size_type = typename LFSV::Traits::SizeType;
308
309 // Get geometry
310 auto geo = eg.geometry();
311
312 // Determine quadrature order
313 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
314 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
315 const int jac_order = geo.type().isSimplex() ? 0 : 1;
316 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
317
318 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
319
320 // loop over quadrature points
321 for (const auto& ip : quadratureRule(geo,qorder))
322 {
323 auto local = ip.position();
324 auto mu = prm.mu(eg,local);
325 auto rho = prm.rho(eg,local);
326
327 // compute u (if Navier term enabled)
328 std::vector<Range_V> phi_v(lfsv_v.size());
329 Range_V val_u(0.0);
330 if(navier) {
331 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
332 for (size_type i=0; i<lfsu_v.size(); i++)
333 val_u.axpy(x(lfsu_v,i),phi_v[i]);
334 }
335
336 // values of pressure shape functions
337 std::vector<Range_P> phi_p(lfsv_p.size());
338 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
339
340 // evaluate jacobian of velocity shape functions on reference element
341 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
342 VectorBasisSwitch_V::jacobian
343 (FESwitch_V::basis(lfsv_v.finiteElement()), geo, local, jac_phi_v);
344
345 assert(lfsu_v.size() == lfsv_v.size());
346 // compute divergence of velocity shape functions
347 std::vector<RF> div_phi_v(lfsv_v.size(),0.0);
348 for (size_type i=0; i<lfsv_v.size(); i++)
349 for (size_type d=0; d<dim; d++)
350 div_phi_v[i] += jac_phi_v[i][d][d];
351
352 // compute velocity jacobian (if Navier term enabled)
354 if(navier) {
355 for (size_type i=0; i<lfsu_v.size(); i++){
356 jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
357 }
358 }
359
360 auto detj = geo.integrationElement(ip.position());
361 auto weight = ip.weight() * detj;
362
363 for(size_type i=0; i<lfsv_v.size(); i++) {
364
365 for(size_type j=0; j<lfsu_v.size(); j++) {
366 //================================================//
367 // \int (mu*grad_u*grad_v)
368 //================================================//
369 RF dvdu(0.0);
370 contraction(jac_phi_v[j],jac_phi_v[i],dvdu);
371 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * dvdu * weight);
372
373 //================================================//
374 // \int \rho ((u\cdot\nabla ) u )\cdot v
375 //================================================//
376 if(navier) {
377 // compute (grad u) phi_v_j (matrix-vector product)
378 Range_V nabla_u_phi_v_j(0.0);
379 jac_u.mv(phi_v[j],nabla_u_phi_v_j);
380 // compute (grad phi_v_j) u (matrix-vector product)
381 Range_V nabla_phi_v_j_u(0.0);
382 jac_phi_v[j].mv(val_u,nabla_phi_v_j_u);
383 mat.accumulate(lfsv_v,i,lfsu_v,j, rho * ((nabla_u_phi_v_j*phi_v[i]) + (nabla_phi_v_j_u*phi_v[i])) * weight);
384 } // end navier
385
386 } // end j
387
388 for(size_type j=0; j<lfsv_p.size(); j++) {
389 //================================================//
390 // - p * div v
391 // - q * div u
392 //================================================//
393 mat.accumulate(lfsv_v,i,lfsu_p,j, -phi_p[j] * div_phi_v[i] * weight);
394 mat.accumulate(lfsv_p,j,lfsu_v,i, -phi_p[j] * div_phi_v[i] * incomp_scaling * weight);
395 }
396 } // end i
397
398 } // end loop quadrature points
399 } // end jacobian_volume
400
401 // skeleton term, each face is only visited ONCE
402 template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
403 void alpha_skeleton (const IG& ig,
404 const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
405 const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n,
406 R& r_s, R& r_n) const
407 {
408 // dimensions
409 const unsigned int dim = IG::Entity::dimension;
410 const unsigned int dimw = IG::coorddimension;
411
412 // subspaces
413 using namespace Indices;
414 using LFSV_V = TypeTree::Child<LFSV,_0>;
415 const auto& lfsv_v_s = child(lfsv_s,_0);
416 const auto& lfsu_v_s = child(lfsu_s,_0);
417 const auto& lfsv_v_n = child(lfsv_n,_0);
418 const auto& lfsu_v_n = child(lfsu_n,_0);
419
420 using LFSV_P = TypeTree::Child<LFSV,_1>;
421 const auto& lfsv_p_s = child(lfsv_s,_1);
422 const auto& lfsu_p_s = child(lfsu_s,_1);
423 const auto& lfsv_p_n = child(lfsv_n,_1);
424 const auto& lfsu_p_n = child(lfsu_n,_1);
425
426 // domain and range field type
427 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
428 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
429 using VectorBasisSwitch_V = VectorBasisInterfaceSwitch<typename FESwitch_V::Basis >;
430 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
431 using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
432 using DF = typename BasisSwitch_V::DomainField;
433 using RF = typename BasisSwitch_V::RangeField;
434 using Range_V = typename BasisSwitch_V::Range;
435 using Range_P = typename BasisSwitch_P::Range;
436 using size_type = typename LFSV::Traits::SizeType;
437
438 // References to inside and outside cells
439 const auto& cell_inside = ig.inside();
440 const auto& cell_outside = ig.outside();
441
442 // Get geometries
443 auto geo = ig.geometry();
444 auto geo_inside = cell_inside.geometry();
445 auto geo_outside = cell_outside.geometry();
446
447 // Get geometry of intersection in local coordinates of cell_inside and cell_outside
448 auto geo_in_inside = ig.geometryInInside();
449 auto geo_in_outside = ig.geometryInOutside();
450
451 // Determine quadrature order
452 const int v_order = FESwitch_V::basis(lfsv_v_s.finiteElement()).order();
453 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-2);
454 const int qorder = 2*v_order + det_jac_order + superintegration_order;
455
456 const int epsilon = prm.epsilonIPSymmetryFactor();
457 const RF incomp_scaling = prm.incompressibilityScaling(current_dt);
458
459 auto penalty_factor = prm.getFaceIP(geo,geo_inside,geo_outside);
460
461 // loop over quadrature points and integrate normal flux
462 for (const auto& ip : quadratureRule(geo,qorder))
463 {
464
465 // position of quadrature point in local coordinates of element
466 auto local_s = geo_in_inside.global(ip.position());
467 auto local_n = geo_in_outside.global(ip.position());
468
469 // values of velocity shape functions
470 std::vector<Range_V> phi_v_s(lfsv_v_s.size());
471 std::vector<Range_V> phi_v_n(lfsv_v_n.size());
472 FESwitch_V::basis(lfsv_v_s.finiteElement()).evaluateFunction(local_s,phi_v_s);
473 FESwitch_V::basis(lfsv_v_n.finiteElement()).evaluateFunction(local_n,phi_v_n);
474
475 // values of pressure shape functions
476 std::vector<Range_P> phi_p_s(lfsv_p_s.size());
477 std::vector<Range_P> phi_p_n(lfsv_p_n.size());
478 FESwitch_P::basis(lfsv_p_s.finiteElement()).evaluateFunction(local_s,phi_p_s);
479 FESwitch_P::basis(lfsv_p_n.finiteElement()).evaluateFunction(local_n,phi_p_n);
480
481 // compute pressure value
482 Range_P val_p_s(0.0);
483 Range_P val_p_n(0.0);
484 for (size_type i=0; i<lfsu_p_s.size(); i++)
485 val_p_s.axpy(x_s(lfsu_p_s,i),phi_p_s[i]);
486 for (size_type i=0; i<lfsu_p_n.size(); i++)
487 val_p_n.axpy(x_n(lfsu_p_n,i),phi_p_n[i]);
488
489 // evaluate jacobian of velocity shape functions on reference element
490 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_s(lfsu_v_s.size());
491 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_n(lfsu_v_n.size());
492 VectorBasisSwitch_V::jacobian
493 (FESwitch_V::basis(lfsv_v_s.finiteElement()), geo_inside, local_s, jac_phi_v_s);
494 VectorBasisSwitch_V::jacobian
495 (FESwitch_V::basis(lfsv_v_n.finiteElement()), geo_outside, local_n, jac_phi_v_n);
496
497 // compute velocity value, jacobian, and divergence
498 Range_V val_u_s(0.0);
499 Range_V val_u_n(0.0);
502 for (size_type i=0; i<lfsu_v_s.size(); i++){
503 val_u_s.axpy(x_s(lfsu_v_s,i),phi_v_s[i]);
504 jac_u_s.axpy(x_s(lfsu_v_s,i),jac_phi_v_s[i]);
505 }
506 for (size_type i=0; i<lfsu_v_n.size(); i++){
507 val_u_n.axpy(x_n(lfsu_v_n,i),phi_v_n[i]);
508 jac_u_n.axpy(x_n(lfsu_v_n,i),jac_phi_v_n[i]);
509 }
510
511 auto normal = ig.unitOuterNormal(ip.position());
512 auto weight = ip.weight()*geo.integrationElement(ip.position());
513 auto mu = prm.mu(ig,ip.position());
514
515 auto factor = mu * weight;
516
517 // compute jump in velocity
518 auto jump = val_u_s - val_u_n;
519
520 // compute mean in pressure
521 auto mean_p = 0.5*(val_p_s + val_p_n);
522
523 // compute flux of velocity jacobian
524 Dune::FieldVector<DF,dimw> flux_jac_u(0.0);
525 add_compute_flux(jac_u_s,normal,flux_jac_u);
526 add_compute_flux(jac_u_n,normal,flux_jac_u);
527 flux_jac_u *= 0.5;
528
529 // loop over test functions, same element
530 for (size_t i=0; i<lfsv_v_s.size(); i++) {
531 //================================================//
532 // diffusion term
533 //================================================//
534 r_s.accumulate(lfsv_v_s, i, -(flux_jac_u * phi_v_s[i]) * factor);
535
536 //================================================//
537 // (non-)symmetric IP term
538 //================================================//
539 Dune::FieldVector<DF,dimw> flux_jac_phi(0.0);
540 add_compute_flux(jac_phi_v_s[i],normal,flux_jac_phi);
541 r_s.accumulate(lfsv_v_s, i, epsilon * 0.5 * (flux_jac_phi * jump) * factor);
542
543 //================================================//
544 // standard IP term integral
545 //================================================//
546 r_s.accumulate(lfsv_v_s,i, penalty_factor * (jump*phi_v_s[i]) * factor);
547
548 //================================================//
549 // pressure-velocity-coupling in momentum equation
550 //================================================//
551 r_s.accumulate(lfsv_v_s,i, mean_p * (phi_v_s[i]*normal) * weight);
552 }
553
554 // loop over test functions, neighbour element
555 for (size_t i=0; i<lfsv_v_n.size(); i++) {
556 //================================================//
557 // diffusion term
558 //================================================//
559 r_n.accumulate(lfsv_v_n, i, (flux_jac_u * phi_v_n[i]) * factor);
560
561 //================================================//
562 // (non-)symmetric IP term
563 //================================================//
564 Dune::FieldVector<DF,dimw> flux_jac_phi(0.0);
565 add_compute_flux(jac_phi_v_n[i],normal,flux_jac_phi);
566 r_n.accumulate(lfsv_v_n, i, epsilon * 0.5 * (flux_jac_phi * jump) * factor);
567
568 //================================================//
569 // standard IP term integral
570 //================================================//
571 r_n.accumulate(lfsv_v_n,i, -penalty_factor * (jump*phi_v_n[i]) * factor);
572
573 //================================================//
574 // pressure-velocity-coupling in momentum equation
575 //================================================//
576 r_n.accumulate(lfsv_v_n,i, -mean_p * (phi_v_n[i]*normal) * weight);
577 }
578
579 //================================================//
580 // incompressibility constraint
581 //================================================//
582 for (size_t i=0; i<lfsv_p_s.size(); i++)
583 r_s.accumulate(lfsv_p_s,i, 0.5*phi_p_s[i] * (jump*normal) * incomp_scaling * weight);
584 for (size_t i=0; i<lfsv_p_n.size(); i++)
585 r_n.accumulate(lfsv_p_n,i, 0.5*phi_p_n[i] * (jump*normal) * incomp_scaling * weight);
586
587 } // end loop quadrature points
588 } // end alpha_skeleton
589
590 // jacobian of skeleton term, each face is only visited ONCE
591 template<typename IG, typename LFSU, typename X, typename LFSV,
592 typename LocalMatrix>
593 void jacobian_skeleton (const IG& ig,
594 const LFSU& lfsu_s, const X&, const LFSV& lfsv_s,
595 const LFSU& lfsu_n, const X&, const LFSV& lfsv_n,
596 LocalMatrix& mat_ss, LocalMatrix& mat_sn,
597 LocalMatrix& mat_ns, LocalMatrix& mat_nn) const
598 {
599 // dimensions
600 const unsigned int dim = IG::Entity::dimension;
601 const unsigned int dimw = IG::coorddimension;
602
603 // subspaces
604 using namespace Indices;
605 using LFSV_V = TypeTree::Child<LFSV,_0>;
606 const auto& lfsv_v_s = child(lfsv_s,_0);
607 const auto& lfsu_v_s = child(lfsu_s,_0);
608 const auto& lfsv_v_n = child(lfsv_n,_0);
609 const auto& lfsu_v_n = child(lfsu_n,_0);
610
611 using LFSV_P = TypeTree::Child<LFSV,_1>;
612 const auto& lfsv_p_s = child(lfsv_s,_1);
613 const auto& lfsu_p_s = child(lfsu_s,_1);
614 const auto& lfsv_p_n = child(lfsv_n,_1);
615 const auto& lfsu_p_n = child(lfsu_n,_1);
616
617 // domain and range field type
618 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
619 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
620 using VectorBasisSwitch_V = VectorBasisInterfaceSwitch<typename FESwitch_V::Basis >;
621 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
622 using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
623 using DF = typename BasisSwitch_V::DomainField;
624 using RF = typename BasisSwitch_V::RangeField;
625 using Range_V = typename BasisSwitch_V::Range;
626 using Range_P = typename BasisSwitch_P::Range;
627 using size_type = typename LFSV::Traits::SizeType;
628
629 // References to inside and outside cells
630 auto const& cell_inside = ig.inside();
631 auto const& cell_outside = ig.outside();
632
633 // Get geometries
634 auto geo = ig.geometry();
635 auto geo_inside = cell_inside.geometry();
636 auto geo_outside = cell_outside.geometry();
637
638 // Get geometry of intersection in local coordinates of cell_inside and cell_outside
639 auto geo_in_inside = ig.geometryInInside();
640 auto geo_in_outside = ig.geometryInOutside();
641
642 // Determine quadrature order
643 const int v_order = FESwitch_V::basis(lfsv_v_s.finiteElement()).order();
644 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-2);
645 const int qorder = 2*v_order + det_jac_order + superintegration_order;
646
647 auto epsilon = prm.epsilonIPSymmetryFactor();
648 auto incomp_scaling = prm.incompressibilityScaling(current_dt);
649
650 auto penalty_factor = prm.getFaceIP(geo,geo_inside,geo_outside);
651
652 // loop over quadrature points and integrate normal flux
653 for (const auto& ip : quadratureRule(geo,qorder))
654 {
655
656 // position of quadrature point in local coordinates of element
657 auto local_s = geo_in_inside.global(ip.position());
658 auto local_n = geo_in_outside.global(ip.position());
659
660 // values of velocity shape functions
661 std::vector<Range_V> phi_v_s(lfsv_v_s.size());
662 std::vector<Range_V> phi_v_n(lfsv_v_n.size());
663 FESwitch_V::basis(lfsv_v_s.finiteElement()).evaluateFunction(local_s,phi_v_s);
664 FESwitch_V::basis(lfsv_v_n.finiteElement()).evaluateFunction(local_n,phi_v_n);
665
666 // values of pressure shape functions
667 std::vector<Range_P> phi_p_s(lfsv_p_s.size());
668 std::vector<Range_P> phi_p_n(lfsv_p_n.size());
669 FESwitch_P::basis(lfsv_p_s.finiteElement()).evaluateFunction(local_s,phi_p_s);
670 FESwitch_P::basis(lfsv_p_n.finiteElement()).evaluateFunction(local_n,phi_p_n);
671
672 // evaluate jacobian of velocity shape functions on reference element
673 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_s(lfsu_v_s.size());
674 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v_n(lfsu_v_n.size());
675 VectorBasisSwitch_V::jacobian
676 (FESwitch_V::basis(lfsv_v_s.finiteElement()), geo_inside, local_s, jac_phi_v_s);
677 VectorBasisSwitch_V::jacobian
678 (FESwitch_V::basis(lfsv_v_n.finiteElement()), geo_outside, local_n, jac_phi_v_n);
679
680 auto normal = ig.unitOuterNormal(ip.position());
681 auto weight = ip.weight()*geo.integrationElement(ip.position());
682 auto mu = prm.mu(ig,ip.position());
683
684 auto factor = mu * weight;
685
686 //============================================
687 // loop over test functions, same element
688 //============================================
689 for(size_type i=0; i<lfsv_v_s.size(); i++) {
690
691 // compute flux
692 Dune::FieldVector<DF,dimw> flux_jac_phi_i(0.0);
693 add_compute_flux(jac_phi_v_s[i],normal,flux_jac_phi_i);
694
695 //============================================
696 // diffusion
697 // (non-)symmetric IP-Term
698 // standard IP integral
699 //============================================
700 for(size_type j=0; j<lfsu_v_s.size(); j++) {
701 Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
702 add_compute_flux(jac_phi_v_s[j],normal,flux_jac_phi_j);
703
704 mat_ss.accumulate(lfsv_v_s,i,lfsu_v_s,j, -0.5 * (flux_jac_phi_j*phi_v_s[i]) * factor);
705 mat_ss.accumulate(lfsv_v_s,i,lfsu_v_s,j, epsilon * 0.5 * (flux_jac_phi_i*phi_v_s[j]) * factor);
706 mat_ss.accumulate(lfsv_v_s,i,lfsu_v_s,j, penalty_factor * (phi_v_s[j]*phi_v_s[i]) * factor);
707 }
708
709 for(size_type j=0; j<lfsu_v_n.size(); j++) {
710 Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
711 add_compute_flux(jac_phi_v_n[j],normal,flux_jac_phi_j);
712
713 mat_sn.accumulate(lfsv_v_s,i,lfsu_v_n,j, -0.5 * (flux_jac_phi_j*phi_v_s[i]) * factor);
714 mat_sn.accumulate(lfsv_v_s,i,lfsu_v_n,j, -epsilon * 0.5 * (flux_jac_phi_i*phi_v_n[j]) * factor);
715 mat_sn.accumulate(lfsv_v_s,i,lfsu_v_n,j, -penalty_factor * (phi_v_n[j]*phi_v_s[i]) * factor);
716 }
717
718 //============================================
719 // pressure-velocity coupling in momentum equation
720 //============================================
721 for(size_type j=0; j<lfsu_p_s.size(); j++) {
722 mat_ss.accumulate(lfsv_v_s,i,lfsu_p_s,j, 0.5*phi_p_s[j] * (phi_v_s[i]*normal) * weight);
723 }
724
725 for(size_type j=0; j<lfsu_p_n.size(); j++) {
726 mat_sn.accumulate(lfsv_v_s,i,lfsu_p_n,j, 0.5*phi_p_n[j] * (phi_v_s[i]*normal) * weight);
727 }
728 } // end i (same)
729
730 //============================================
731 // loop over test functions, neighbour element
732 //============================================
733 for(size_type i=0; i<lfsv_v_n.size(); i++) {
734
735 // compute flux
736 Dune::FieldVector<DF,dimw> flux_jac_phi_i(0.0);
737 add_compute_flux(jac_phi_v_n[i],normal,flux_jac_phi_i);
738
739 //============================================
740 // diffusion
741 // (non-)symmetric IP-Term
742 // standard IP integral
743 //============================================
744 for(size_type j=0; j<lfsu_v_s.size(); j++) {
745 Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
746 add_compute_flux(jac_phi_v_s[j],normal,flux_jac_phi_j);
747
748 mat_ns.accumulate(lfsv_v_n,i,lfsu_v_s,j, 0.5 * (flux_jac_phi_j*phi_v_n[i]) * factor);
749 mat_ns.accumulate(lfsv_v_n,i,lfsu_v_s,j, epsilon * 0.5 * (flux_jac_phi_i*phi_v_s[j]) * factor);
750 mat_ns.accumulate(lfsv_v_n,i,lfsu_v_s,j, -penalty_factor * (phi_v_s[j]*phi_v_n[i]) * factor);
751 }
752
753 for(size_type j=0; j<lfsu_v_n.size(); j++) {
754 Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
755 add_compute_flux(jac_phi_v_n[j],normal,flux_jac_phi_j);
756
757 mat_nn.accumulate(lfsv_v_n,i,lfsu_v_n,j, 0.5 * (flux_jac_phi_j*phi_v_n[i]) * factor);
758 mat_nn.accumulate(lfsv_v_n,i,lfsu_v_n,j, -epsilon * 0.5 * (flux_jac_phi_i*phi_v_n[j]) * factor);
759 mat_nn.accumulate(lfsv_v_n,i,lfsu_v_n,j, penalty_factor * (phi_v_n[j]*phi_v_n[i]) * factor);
760 }
761
762 //============================================
763 // pressure-velocity coupling in momentum equation
764 //============================================
765 for(size_type j=0; j<lfsu_p_s.size(); j++) {
766 mat_ns.accumulate(lfsv_v_n,i,lfsu_p_s,j, -0.5*phi_p_s[j] * (phi_v_n[i]*normal) * weight);
767 }
768
769 for(size_type j=0; j<lfsu_p_n.size(); j++) {
770 mat_nn.accumulate(lfsv_v_n,i,lfsu_p_n,j, -0.5*phi_p_n[j] * (phi_v_n[i]*normal) * weight);
771 }
772 } // end i (neighbour)
773
774 //================================================//
775 // \int <q> [u] n
776 //================================================//
777 for(size_type i=0; i<lfsv_p_s.size(); i++) {
778 for(size_type j=0; j<lfsu_v_s.size(); j++)
779 mat_ss.accumulate(lfsv_p_s,i,lfsu_v_s,j, 0.5*phi_p_s[i] * (phi_v_s[j]*normal) * incomp_scaling * weight);
780
781 for(size_type j=0; j<lfsu_v_n.size(); j++)
782 mat_sn.accumulate(lfsv_p_s,i,lfsu_v_n,j, -0.5*phi_p_s[i] * (phi_v_n[j]*normal) * incomp_scaling * weight);
783 }
784
785 for(size_type i=0; i<lfsv_p_n.size(); i++) {
786 for(size_type j=0; j<lfsu_v_s.size(); j++)
787 mat_ns.accumulate(lfsv_p_n,i,lfsu_v_s,j, 0.5*phi_p_n[i] * (phi_v_s[j]*normal) * incomp_scaling * weight);
788
789 for(size_type j=0; j<lfsu_v_n.size(); j++)
790 mat_nn.accumulate(lfsv_p_n,i,lfsu_v_n,j, -0.5*phi_p_n[i] * (phi_v_n[j]*normal) * incomp_scaling * weight);
791 }
792
793 } // end loop quadrature points
794 } // end jacobian_skeleton
795
796 // boundary term
797 template<typename IG, typename LFSU, typename X, typename LFSV, typename R>
798 void alpha_boundary (const IG& ig,
799 const LFSU& lfsu, const X& x, const LFSV& lfsv,
800 R& r) const
801 {
802 // dimensions
803 const unsigned int dim = IG::Entity::dimension;
804 const unsigned int dimw = IG::coorddimension;
805
806 // subspaces
807 using namespace Indices;
808 using LFSV_V = TypeTree::Child<LFSV,_0>;
809 const auto& lfsv_v = child(lfsv,_0);
810 const auto& lfsu_v = child(lfsu,_0);
811
812 using LFSV_P = TypeTree::Child<LFSV,_1>;
813 const auto& lfsv_p = child(lfsv,_1);
814 const auto& lfsu_p = child(lfsu,_1);
815
816 // domain and range field type
817 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
818 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
819 using VectorBasisSwitch_V = VectorBasisInterfaceSwitch<typename FESwitch_V::Basis >;
820 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
821 using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
822 using DF = typename BasisSwitch_V::DomainField;
823 using RF = typename BasisSwitch_V::RangeField;
824 using Range_V = typename BasisSwitch_V::Range;
825 using Range_P = typename BasisSwitch_P::Range;
826 using size_type = typename LFSV::Traits::SizeType;
827
828 // References to inside cell
829 const auto& cell_inside = ig.inside();
830
831 // Get geometries
832 auto geo = ig.geometry();
833 auto geo_inside = cell_inside.geometry();
834
835 // Get geometry of intersection in local coordinates of cell_inside
836 auto geo_in_inside = ig.geometryInInside();
837
838 // Determine quadrature order
839 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
840 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
841 const int qorder = 2*v_order + det_jac_order + superintegration_order;
842
843 auto epsilon = prm.epsilonIPSymmetryFactor();
844 auto incomp_scaling = prm.incompressibilityScaling(current_dt);
845
846 auto penalty_factor = prm.getFaceIP(geo,geo_inside);
847
848 // loop over quadrature points and integrate normal flux
849 for (const auto& ip : quadratureRule(geo,qorder))
850 {
851 // position of quadrature point in local coordinates of element
852 auto local = geo_in_inside.global(ip.position());
853
854 // values of velocity shape functions
855 std::vector<Range_V> phi_v(lfsv_v.size());
856 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
857
858 // values of pressure shape functions
859 std::vector<Range_P> phi_p(lfsv_p.size());
860 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
861
862 // evaluate jacobian of basis functions on reference element
863 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
864 VectorBasisSwitch_V::jacobian
865 (FESwitch_V::basis(lfsv_v.finiteElement()), geo_inside, local, jac_phi_v);
866
867 // compute pressure value
868 Range_P val_p(0.0);
869 for (size_type i=0; i<lfsu_p.size(); i++)
870 val_p.axpy(x(lfsu_p,i),phi_p[i]);
871
872 // compute u and velocity jacobian
873 Range_V val_u(0.0);
875 for (size_type i=0; i<lfsu_v.size(); i++){
876 val_u.axpy(x(lfsu_v,i),phi_v[i]);
877 jac_u.axpy(x(lfsu_v,i),jac_phi_v[i]);
878 }
879
880 auto normal = ig.unitOuterNormal(ip.position());
881 auto weight = ip.weight()*geo.integrationElement(ip.position());
882 auto mu = prm.mu(ig,ip.position());
883
884 // evaluate boundary condition type
885 auto bctype(prm.bctype(ig,ip.position()));
886
887 if (bctype == BC::VelocityDirichlet) {
888 // compute jump relative to Dirichlet value
889 auto u0(prm.g(cell_inside,local));
890 auto jump = val_u - u0;
891
892 // compute flux of velocity jacobian
893 Dune::FieldVector<DF,dimw> flux_jac_u(0.0);
894 add_compute_flux(jac_u,normal,flux_jac_u);
895
896 for (size_t i=0; i<lfsv_v.size(); i++) {
897 //================================================//
898 // diffusion term
899 //================================================//
900 r.accumulate(lfsv_v,i, -mu * (flux_jac_u * phi_v[i]) * weight);
901
902 //================================================//
903 // (non-)symmetric IP term
904 //================================================//
905 Dune::FieldVector<DF,dimw> flux_jac_phi(0.0);
906 add_compute_flux(jac_phi_v[i],normal,flux_jac_phi);
907 r.accumulate(lfsv_v,i, mu * epsilon * (flux_jac_phi*jump) * weight);
908
909 //================================================//
910 // standard IP term integral
911 //================================================//
912 r.accumulate(lfsv_v,i, mu * (jump*phi_v[i]) * penalty_factor * weight);
913
914 //================================================//
915 // pressure-velocity-coupling in momentum equation
916 //================================================//
917 r.accumulate(lfsv_v,i, val_p * (phi_v[i]*normal) * weight);
918 } // end i
919
920 //================================================//
921 // incompressibility constraint
922 //================================================//
923 for(size_type i=0; i<lfsv_p.size(); i++) {
924 r.accumulate(lfsv_p,i, phi_p[i] * (jump*normal) * incomp_scaling * weight);
925 }
926 } // Velocity Dirichlet
927
928 if (bctype == BC::StressNeumann) {
929 auto stress(prm.j(ig,ip.position(),normal));
930
931 for(size_type i=0; i<lfsv_v.size(); i++) {
932 r.accumulate(lfsv_v,i, (stress*phi_v[i]) * weight);
933 }
934 } // Pressure Dirichlet
935
936 } // end loop quadrature points
937 } // end alpha_boundary
938
939 // jacobian of boundary term
940 template<typename IG, typename LFSU, typename X, typename LFSV,
941 typename LocalMatrix>
942 void jacobian_boundary (const IG& ig,
943 const LFSU& lfsu, const X& x, const LFSV& lfsv,
944 LocalMatrix& mat) const
945 {
946 // dimensions
947 const unsigned int dim = IG::Entity::dimension;
948 const unsigned int dimw = IG::coorddimension;
949
950 // subspaces
951 using namespace Indices;
952 using LFSV_V = TypeTree::Child<LFSV,_0>;
953 const auto& lfsv_v = child(lfsv,_0);
954 const auto& lfsu_v = child(lfsu,_0);
955
956 using LFSV_P = TypeTree::Child<LFSV,_1>;
957 const auto& lfsv_p = child(lfsv,_1);
958 const auto& lfsu_p = child(lfsu,_1);
959
960 // domain and range field type
961 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
962 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
963 using VectorBasisSwitch_V = VectorBasisInterfaceSwitch<typename FESwitch_V::Basis >;
964 using FESwitch_P = FiniteElementInterfaceSwitch<typename LFSV_P::Traits::FiniteElementType >;
965 using BasisSwitch_P = BasisInterfaceSwitch<typename FESwitch_P::Basis >;
966 using DF = typename BasisSwitch_V::DomainField;
967 using RF = typename BasisSwitch_V::RangeField;
968 using Range_V = typename BasisSwitch_V::Range;
969 using Range_P = typename BasisSwitch_P::Range;
970 using size_type = typename LFSV::Traits::SizeType;
971
972 // References to inside cell
973 const auto& cell_inside = ig.inside();
974
975 // Get geometries
976 auto geo = ig.geometry();
977 auto geo_inside = cell_inside.geometry();
978
979 // Get geometry of intersection in local coordinates of cell_inside
980 auto geo_in_inside = ig.geometryInInside();
981
982 // Determine quadrature order
983 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
984 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
985 const int qorder = 2*v_order + det_jac_order + superintegration_order;
986
987 auto epsilon = prm.epsilonIPSymmetryFactor();
988 auto incomp_scaling = prm.incompressibilityScaling(current_dt);
989
990 auto penalty_factor = prm.getFaceIP(geo,geo_inside);
991
992 // loop over quadrature points and integrate normal flux
993 for (const auto& ip : quadratureRule(geo,qorder))
994 {
995 // position of quadrature point in local coordinates of element
996 auto local = geo_in_inside.global(ip.position());
997
998 // values of velocity shape functions
999 std::vector<Range_V> phi_v(lfsv_v.size());
1000 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1001
1002 // values of pressure shape functions
1003 std::vector<Range_P> phi_p(lfsv_p.size());
1004 FESwitch_P::basis(lfsv_p.finiteElement()).evaluateFunction(local,phi_p);
1005
1006 // evaluate jacobian of basis functions on reference element
1007 std::vector<Dune::FieldMatrix<RF,dim,dim> > jac_phi_v(lfsu_v.size());
1008 VectorBasisSwitch_V::jacobian
1009 (FESwitch_V::basis(lfsv_v.finiteElement()), geo_inside, local, jac_phi_v);
1010
1011 auto normal = ig.unitOuterNormal(ip.position());
1012 auto weight = ip.weight()*geo.integrationElement(ip.position());
1013 auto mu = prm.mu(ig,ip.position());
1014
1015 // evaluate boundary condition type
1016 auto bctype(prm.bctype(ig,ip.position()));
1017
1018 if (bctype == BC::VelocityDirichlet) {
1019
1020 for(size_type i=0; i<lfsv_v.size(); i++) {
1021 // compute flux
1022 Dune::FieldVector<DF,dimw> flux_jac_phi_i(0.0);
1023 add_compute_flux(jac_phi_v[i],normal,flux_jac_phi_i);
1024
1025 for(size_type j=0; j<lfsu_v.size(); j++) {
1026 //================================================//
1027 // diffusion term
1028 // (non-)symmetric IP term
1029 //================================================//
1030 Dune::FieldVector<DF,dimw> flux_jac_phi_j(0.0);
1031 add_compute_flux(jac_phi_v[j],normal,flux_jac_phi_j);
1032
1033 mat.accumulate(lfsv_v,i,lfsu_v,j, -mu * (flux_jac_phi_j*phi_v[i]) * weight);
1034 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * epsilon * (flux_jac_phi_i*phi_v[j]) *weight);
1035
1036 //================================================//
1037 // standard IP term integral
1038 //================================================//
1039 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (phi_v[j]*phi_v[i]) * penalty_factor * weight);
1040 }
1041
1042 //================================================//
1043 // pressure-velocity-coupling in momentum equation
1044 //================================================//
1045 for(size_type j=0; j<lfsu_p.size(); j++) {
1046 mat.accumulate(lfsv_v,i,lfsu_p,j, phi_p[j] * (phi_v[i]*normal) * weight);
1047 }
1048 } // end i
1049
1050 //================================================//
1051 // incompressibility constraint
1052 //================================================//
1053 for(size_type i=0; i<lfsv_p.size(); i++) {
1054 for(size_type j=0; j<lfsu_v.size(); j++) {
1055 mat.accumulate(lfsv_p,i,lfsu_v,j, phi_p[i] * (phi_v[j]*normal) * incomp_scaling * weight);
1056 }
1057 }
1058
1059 } // Velocity Dirichlet
1060
1061 } // end loop quadrature points
1062 } // end jacobian_boundary
1063
1064 // volume integral depending only on test functions,
1065 // contains f on the right hand side
1066 template<typename EG, typename LFSV, typename R>
1067 void lambda_volume (const EG& eg, const LFSV& lfsv, R& r) const
1068 {
1069 const unsigned int dim = EG::Geometry::mydimension;
1070
1071 // subspaces
1072 using namespace Indices;
1073 using LFSV_V = TypeTree::Child<LFSV,_0>;
1074 const auto& lfsv_v = child(lfsv,_0);
1075
1076 // domain and range field type
1077 using FESwitch_V = FiniteElementInterfaceSwitch<typename LFSV_V::Traits::FiniteElementType >;
1078 using BasisSwitch_V = BasisInterfaceSwitch<typename FESwitch_V::Basis >;
1079 using Range_V = typename BasisSwitch_V::Range;
1080 using size_type = typename LFSV::Traits::SizeType;
1081
1082 // Get cell
1083 const auto& cell = eg.entity();
1084
1085 // Get geometries
1086 auto geo = eg.geometry();
1087
1088 // Determine quadrature order
1089 const int v_order = FESwitch_V::basis(lfsv_v.finiteElement()).order();
1090 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
1091 const int qorder = 2*v_order + det_jac_order + superintegration_order;
1092
1093 // loop over quadrature points
1094 for (const auto& ip : quadratureRule(geo,qorder))
1095 {
1096 auto local = ip.position();
1097 //const Dune::FieldVector<DF,dimw> global = eg.geometry().global(local);
1098
1099 // values of velocity shape functions
1100 std::vector<Range_V> phi_v(lfsv_v.size());
1101 FESwitch_V::basis(lfsv_v.finiteElement()).evaluateFunction(local,phi_v);
1102
1103 auto weight = ip.weight() * geo.integrationElement(ip.position());
1104
1105 // evaluate source term
1106 auto fval(prm.f(cell,local));
1107
1108 //================================================//
1109 // \int (f*v)
1110 //================================================//
1111 for(size_type i=0; i<lfsv_v.size(); i++)
1112 r.accumulate(lfsv_v,i, -(fval*phi_v[i]) * weight);
1113
1114 } // end loop quadrature points
1115 } // end lambda_volume
1116
1117 private :
1118
1119 template<class M, class RF>
1120 static void contraction(const M & a, const M & b, RF & v)
1121 {
1122 v = 0;
1123 const int an = a.N();
1124 const int am = a.M();
1125 for(int r=0; r<an; ++r)
1126 for(int c=0; c<am; ++c)
1127 v += a[r][c] * b[r][c];
1128 }
1129
1130 template<class DU, class R>
1131 static void add_compute_flux(const DU & du, const R & n, R & result)
1132 {
1133 const int N = du.N();
1134 const int M = du.M();
1135 for(int r=0; r<N; ++r)
1136 for(int c=0; c<M; ++c)
1137 result[r] += du[r][c] * n[c];
1138 }
1139
1140 PRM& prm;
1141 const int superintegration_order;
1142 Real current_dt;
1143 }; // end class DGNavierStokesVelVecFEM
1144
1145 } // end namespace PDELab
1146} // end namespace Dune
1147#endif // DUNE_PDELAB_LOCALOPERATOR_DGNAVIERSTOKESVELVECFEM_HH
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Wrapper class for geometries.
Definition: geometry.hh:71
Implementation::JacobianInverseTransposed JacobianInverseTransposed
type of jacobian inverse transposed
Definition: geometry.hh:120
static constexpr int coorddimension
dimension of embedding coordinate system
Definition: geometry.hh:97
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
Return inverse of transposed of Jacobian.
Definition: geometry.hh:328
A local operator for solving the Navier-Stokes equations using a DG discretization and a vector-value...
Definition: dgnavierstokesvelvecfem.hh:105
DGNavierStokesVelVecFEM(PRM &_prm, int _superintegration_order=0)
Constructor.
Definition: dgnavierstokesvelvecfem.hh:142
sparsity pattern generator
Definition: pattern.hh:30
sparsity pattern generator
Definition: pattern.hh:14
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
void setTime(R t_)
set time for subsequent evaluation
Definition: idefault.hh:104
Default flags for all local operators.
Definition: flags.hh:19
A few common exception classes.
Traits for type conversions and type information.
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
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
STL namespace.
A hierarchical structure of string parameters.
Definition: stokesparameter.hh:32
static void jacobian(const Basis &basis, const Geometry &geometry, const DomainLocal &xl, std::vector< FieldMatrix< RangeField, dimRange, Geometry::coorddimension > > &jac)
Compute global jacobian matrix for vector valued bases.
Definition: dgnavierstokesvelvecfem.hh:73
typename Basis::Traits::DomainType DomainLocal
export vector type of the local coordinates
Definition: dgnavierstokesvelvecfem.hh:65
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)