2 #ifndef DUNE_PDELAB_CONVECTIONDIFFUSIONDG_HH
3 #define DUNE_PDELAB_CONVECTIONDIFFUSIONDG_HH
5 #include<dune/common/exceptions.hh>
6 #include<dune/common/fvector.hh>
7 #include<dune/geometry/referenceelements.hh>
56 template<
typename T,
typename FiniteElementMap>
66 enum { dim = T::Traits::GridViewType::dimension };
68 typedef typename T::Traits::RangeFieldType Real;
92 param(param_), method(method_), weights(weights_),
93 alpha(alpha_), intorderadd(intorderadd_), quadrature_factor(2),
101 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
102 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
105 typedef typename LFSU::Traits::FiniteElementType::
106 Traits::LocalBasisType::Traits::DomainFieldType DF;
107 typedef typename LFSU::Traits::FiniteElementType::
108 Traits::LocalBasisType::Traits::RangeFieldType RF;
109 typedef typename LFSU::Traits::FiniteElementType::
110 Traits::LocalBasisType::Traits::JacobianType JacobianType;
111 typedef typename LFSU::Traits::FiniteElementType::
112 Traits::LocalBasisType::Traits::RangeType RangeType;
113 typedef typename LFSU::Traits::SizeType size_type;
116 const int dim = EG::Entity::dimension;
117 const int order = std::max(lfsu.finiteElement().localBasis().order(),
118 lfsv.finiteElement().localBasis().order());
119 const int intorder = intorderadd + quadrature_factor * order;
122 Dune::GeometryType gt = eg.geometry().type();
123 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
126 typename T::Traits::PermTensorType A;
127 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
128 A = param.A(eg.entity(),localcenter);
131 typename EG::Geometry::JacobianInverseTransposed jac;
134 for (
const auto& ip : rule)
138 std::vector<RangeType> phi(lfsu.size());
139 lfsu.finiteElement().localBasis().evaluateFunction(ip.position(),phi);
140 std::vector<RangeType> psi(lfsv.size());
141 lfsv.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
143 const std::vector<RangeType>& phi = cache[order].evaluateFunction(ip.position(),lfsu.finiteElement().localBasis());
144 const std::vector<RangeType>& psi = cache[order].evaluateFunction(ip.position(),lfsv.finiteElement().localBasis());
149 for (size_type i=0; i<lfsu.size(); i++)
150 u += x(lfsu,i)*phi[i];
154 std::vector<JacobianType> js(lfsu.size());
155 lfsu.finiteElement().localBasis().evaluateJacobian(ip.position(),js);
156 std::vector<JacobianType> js_v(lfsv.size());
157 lfsv.finiteElement().localBasis().evaluateJacobian(ip.position(),js_v);
159 const std::vector<JacobianType>& js = cache[order].evaluateJacobian(ip.position(),lfsu.finiteElement().localBasis());
160 const std::vector<JacobianType>& js_v = cache[order].evaluateJacobian(ip.position(),lfsv.finiteElement().localBasis());
164 jac = eg.geometry().jacobianInverseTransposed(ip.position());
165 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
166 for (size_type i=0; i<lfsu.size(); i++)
167 jac.mv(js[i][0],gradphi[i]);
169 std::vector<Dune::FieldVector<RF,dim> > gradpsi(lfsv.size());
170 for (size_type i=0; i<lfsv.size(); i++)
171 jac.mv(js_v[i][0],gradpsi[i]);
174 Dune::FieldVector<RF,dim> gradu(0.0);
175 for (size_type i=0; i<lfsu.size(); i++)
176 gradu.axpy(x(lfsu,i),gradphi[i]);
179 Dune::FieldVector<RF,dim> Agradu(0.0);
183 typename T::Traits::RangeType b = param.b(eg.entity(),ip.position());
186 typename T::Traits::RangeFieldType c = param.c(eg.entity(),ip.position());
189 RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
190 for (size_type i=0; i<lfsv.size(); i++)
191 r.accumulate(lfsv,i,( Agradu*gradpsi[i] - u*(b*gradpsi[i]) + c*u*psi[i] )*factor);
196 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
201 typedef typename LFSU::Traits::FiniteElementType::
202 Traits::LocalBasisType::Traits::DomainFieldType DF;
203 typedef typename LFSU::Traits::FiniteElementType::
204 Traits::LocalBasisType::Traits::RangeFieldType RF;
205 typedef typename LFSU::Traits::FiniteElementType::
206 Traits::LocalBasisType::Traits::JacobianType JacobianType;
207 typedef typename LFSU::Traits::FiniteElementType::
208 Traits::LocalBasisType::Traits::RangeType RangeType;
209 typedef typename LFSU::Traits::SizeType size_type;
212 const int dim = EG::Entity::dimension;
213 const int order = std::max(lfsu.finiteElement().localBasis().order(),
214 lfsv.finiteElement().localBasis().order());
215 const int intorder = intorderadd + quadrature_factor * order;
218 Dune::GeometryType gt = eg.geometry().type();
219 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
222 typename T::Traits::PermTensorType A;
223 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
224 A = param.A(eg.entity(),localcenter);
227 typename EG::Geometry::JacobianInverseTransposed jac;
230 for (
const auto& ip : rule)
234 std::vector<RangeType> phi(lfsu.size());
235 lfsu.finiteElement().localBasis().evaluateFunction(ip.position(),phi);
237 const std::vector<RangeType>& phi = cache[order].evaluateFunction(ip.position(),lfsu.finiteElement().localBasis());
242 std::vector<JacobianType> js(lfsu.size());
243 lfsu.finiteElement().localBasis().evaluateJacobian(ip.position(),js);
245 const std::vector<JacobianType>& js = cache[order].evaluateJacobian(ip.position(),lfsu.finiteElement().localBasis());
249 jac = eg.geometry().jacobianInverseTransposed(ip.position());
250 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
251 std::vector<Dune::FieldVector<RF,dim> > Agradphi(lfsu.size());
252 for (size_type i=0; i<lfsu.size(); i++)
254 jac.mv(js[i][0],gradphi[i]);
255 A.mv(gradphi[i],Agradphi[i]);
259 typename T::Traits::RangeType b = param.b(eg.entity(),ip.position());
262 typename T::Traits::RangeFieldType c = param.c(eg.entity(),ip.position());
265 RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
266 for (size_type j=0; j<lfsu.size(); j++)
267 for (size_type i=0; i<lfsu.size(); i++)
268 mat.accumulate(lfsu,i,lfsu,j,( Agradphi[j]*gradphi[i] - phi[j]*(b*gradphi[i]) + c*phi[j]*phi[i] )*factor);
274 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
276 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
277 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
278 R& r_s, R& r_n)
const
281 typedef typename LFSV::Traits::FiniteElementType::
282 Traits::LocalBasisType::Traits::DomainFieldType DF;
283 typedef typename LFSV::Traits::FiniteElementType::
284 Traits::LocalBasisType::Traits::RangeFieldType RF;
285 typedef typename LFSV::Traits::FiniteElementType::
286 Traits::LocalBasisType::Traits::RangeType RangeType;
287 typedef typename LFSU::Traits::FiniteElementType::
288 Traits::LocalBasisType::Traits::JacobianType JacobianType;
289 typedef typename LFSV::Traits::SizeType size_type;
292 const int dim = IG::dimension;
293 const int order = std::max(
294 std::max(lfsu_s.finiteElement().localBasis().order(),
295 lfsu_n.finiteElement().localBasis().order()),
296 std::max(lfsv_s.finiteElement().localBasis().order(),
297 lfsv_n.finiteElement().localBasis().order())
299 const int intorder = intorderadd+quadrature_factor*order;
302 auto inside_cell = ig.inside();
303 auto outside_cell = ig.outside();
306 const Dune::FieldVector<DF,dim>&
307 inside_local = Dune::ReferenceElements<DF,dim>::general(inside_cell.type()).position(0,0);
308 const Dune::FieldVector<DF,dim>&
309 outside_local = Dune::ReferenceElements<DF,dim>::general(outside_cell.type()).position(0,0);
310 typename T::Traits::PermTensorType A_s, A_n;
311 A_s = param.A(inside_cell,inside_local);
312 A_n = param.A(outside_cell,outside_local);
318 element_size(inside_cell.geometry(),h_s,hmax_s);
319 element_size(outside_cell.geometry(),h_n,hmax_n);
320 RF h_F = std::min(h_s,h_n);
321 h_F = std::min(inside_cell.geometry().volume(),outside_cell.geometry().volume())/ig.geometry().volume();
324 Dune::GeometryType gtface = ig.geometryInInside().type();
325 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
328 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
331 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
332 Dune::FieldVector<RF,dim> An_F_s;
334 Dune::FieldVector<RF,dim> An_F_n;
340 RF harmonic_average(0.0);
343 RF delta_s = (An_F_s*n_F);
344 RF delta_n = (An_F_n*n_F);
345 omega_s = delta_n/(delta_s+delta_n+1
e-20);
346 omega_n = delta_s/(delta_s+delta_n+1
e-20);
347 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1
e-20);
351 omega_s = omega_n = 0.5;
352 harmonic_average = 1.0;
356 const int order_s = lfsu_s.finiteElement().localBasis().order();
357 const int order_n = lfsu_n.finiteElement().localBasis().order();
358 int degree = std::max( order_s, order_n );
361 RF penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
364 for (
const auto& ip : rule)
367 const Dune::FieldVector<DF,dim> n_F_local = ig.unitOuterNormal(ip.position());
370 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(ip.position());
371 Dune::FieldVector<DF,dim> iplocal_n = ig.geometryInOutside().global(ip.position());
375 std::vector<RangeType> phi_s(lfsu_s.size());
376 lfsu_s.finiteElement().localBasis().evaluateFunction(iplocal_s,phi_s);
377 std::vector<RangeType> phi_n(lfsu_n.size());
378 lfsu_n.finiteElement().localBasis().evaluateFunction(iplocal_n,phi_n);
379 std::vector<RangeType> psi_s(lfsv_s.size());
380 lfsv_s.finiteElement().localBasis().evaluateFunction(iplocal_s,psi_s);
381 std::vector<RangeType> psi_n(lfsv_n.size());
382 lfsv_n.finiteElement().localBasis().evaluateFunction(iplocal_n,psi_n);
384 const std::vector<RangeType>& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
385 const std::vector<RangeType>& phi_n = cache[order_n].evaluateFunction(iplocal_n,lfsu_n.finiteElement().localBasis());
386 const std::vector<RangeType>& psi_s = cache[order_s].evaluateFunction(iplocal_s,lfsv_s.finiteElement().localBasis());
387 const std::vector<RangeType>& psi_n = cache[order_n].evaluateFunction(iplocal_n,lfsv_n.finiteElement().localBasis());
392 for (size_type i=0; i<lfsu_s.size(); i++)
393 u_s += x_s(lfsu_s,i)*phi_s[i];
395 for (size_type i=0; i<lfsu_n.size(); i++)
396 u_n += x_n(lfsu_n,i)*phi_n[i];
400 std::vector<JacobianType> gradphi_s(lfsu_s.size());
401 lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
402 std::vector<JacobianType> gradphi_n(lfsu_n.size());
403 lfsu_n.finiteElement().localBasis().evaluateJacobian(iplocal_n,gradphi_n);
404 std::vector<JacobianType> gradpsi_s(lfsv_s.size());
405 lfsv_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradpsi_s);
406 std::vector<JacobianType> gradpsi_n(lfsv_n.size());
407 lfsv_n.finiteElement().localBasis().evaluateJacobian(iplocal_n,gradpsi_n);
409 const std::vector<JacobianType>& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
410 const std::vector<JacobianType>& gradphi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsu_n.finiteElement().localBasis());
411 const std::vector<JacobianType>& gradpsi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsv_s.finiteElement().localBasis());
412 const std::vector<JacobianType>& gradpsi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsv_n.finiteElement().localBasis());
416 jac = inside_cell.geometry().jacobianInverseTransposed(iplocal_s);
417 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
418 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
419 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_s(lfsv_s.size());
420 for (size_type i=0; i<lfsv_s.size(); i++) jac.mv(gradpsi_s[i][0],tgradpsi_s[i]);
421 jac = outside_cell.geometry().jacobianInverseTransposed(iplocal_n);
422 std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
423 for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
424 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_n(lfsv_n.size());
425 for (size_type i=0; i<lfsv_n.size(); i++) jac.mv(gradpsi_n[i][0],tgradpsi_n[i]);
428 Dune::FieldVector<RF,dim> gradu_s(0.0);
429 for (size_type i=0; i<lfsu_s.size(); i++)
430 gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
431 Dune::FieldVector<RF,dim> gradu_n(0.0);
432 for (size_type i=0; i<lfsu_n.size(); i++)
433 gradu_n.axpy(x_n(lfsu_n,i),tgradphi_n[i]);
436 typename T::Traits::RangeType b = param.b(inside_cell,iplocal_s);
437 RF normalflux = b*n_F_local;
438 RF omegaup_s, omegaup_n;
451 RF factor = ip.weight() * ig.geometry().integrationElement(ip.position());
454 RF term1 = (omegaup_s*u_s + omegaup_n*u_n) * normalflux *factor;
455 for (size_type i=0; i<lfsv_s.size(); i++)
456 r_s.accumulate(lfsu_s,i,term1 * psi_s[i]);
457 for (size_type i=0; i<lfsv_n.size(); i++)
458 r_n.accumulate(lfsu_n,i,-term1 * psi_n[i]);
461 RF term2 = -(omega_s*(An_F_s*gradu_s) + omega_n*(An_F_n*gradu_n)) * factor;
462 for (size_type i=0; i<lfsv_s.size(); i++)
463 r_s.accumulate(lfsv_s,i,term2 * psi_s[i]);
464 for (size_type i=0; i<lfsv_n.size(); i++)
465 r_n.accumulate(lfsv_n,i,-term2 * psi_n[i]);
468 RF term3 = (u_s-u_n) * factor;
469 for (size_type i=0; i<lfsv_s.size(); i++)
470 r_s.accumulate(lfsv_s,i,term3 * theta * omega_s * (An_F_s*tgradpsi_s[i]));
471 for (size_type i=0; i<lfsv_n.size(); i++)
472 r_n.accumulate(lfsv_n,i,term3 * theta * omega_n * (An_F_n*tgradpsi_n[i]));
475 RF term4 = penalty_factor * (u_s-u_n) * factor;
476 for (size_type i=0; i<lfsv_s.size(); i++)
477 r_s.accumulate(lfsv_s,i,term4 * psi_s[i]);
478 for (size_type i=0; i<lfsv_n.size(); i++)
479 r_n.accumulate(lfsv_n,i,-term4 * psi_n[i]);
483 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
485 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
486 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
487 M& mat_ss, M& mat_sn,
488 M& mat_ns, M& mat_nn)
const
491 typedef typename LFSV::Traits::FiniteElementType::
492 Traits::LocalBasisType::Traits::DomainFieldType DF;
493 typedef typename LFSV::Traits::FiniteElementType::
494 Traits::LocalBasisType::Traits::RangeFieldType RF;
495 typedef typename LFSV::Traits::FiniteElementType::
496 Traits::LocalBasisType::Traits::RangeType RangeType;
497 typedef typename LFSU::Traits::FiniteElementType::
498 Traits::LocalBasisType::Traits::JacobianType JacobianType;
499 typedef typename LFSV::Traits::SizeType size_type;
502 const int dim = IG::dimension;
503 const int order = std::max(
504 std::max(lfsu_s.finiteElement().localBasis().order(),
505 lfsu_n.finiteElement().localBasis().order()),
506 std::max(lfsv_s.finiteElement().localBasis().order(),
507 lfsv_n.finiteElement().localBasis().order())
509 const int intorder = intorderadd+quadrature_factor*order;
512 auto inside_cell = ig.inside();
513 auto outside_cell = ig.outside();
516 const Dune::FieldVector<DF,dim>&
517 inside_local = Dune::ReferenceElements<DF,dim>::general(inside_cell.type()).position(0,0);
518 const Dune::FieldVector<DF,dim>&
519 outside_local = Dune::ReferenceElements<DF,dim>::general(outside_cell.type()).position(0,0);
520 typename T::Traits::PermTensorType A_s, A_n;
521 A_s = param.A(inside_cell,inside_local);
522 A_n = param.A(outside_cell,outside_local);
526 DF hmax_s = 0., hmax_n = 0.;
527 element_size(inside_cell.geometry(),h_s,hmax_s);
528 element_size(outside_cell.geometry(),h_n,hmax_n);
529 RF h_F = std::min(h_s,h_n);
530 h_F = std::min(inside_cell.geometry().volume(),outside_cell.geometry().volume())/ig.geometry().volume();
533 Dune::GeometryType gtface = ig.geometryInInside().type();
534 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
537 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
540 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
541 Dune::FieldVector<RF,dim> An_F_s;
543 Dune::FieldVector<RF,dim> An_F_n;
549 RF harmonic_average(0.0);
552 RF delta_s = (An_F_s*n_F);
553 RF delta_n = (An_F_n*n_F);
554 omega_s = delta_n/(delta_s+delta_n+1
e-20);
555 omega_n = delta_s/(delta_s+delta_n+1
e-20);
556 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1
e-20);
560 omega_s = omega_n = 0.5;
561 harmonic_average = 1.0;
565 const int order_s = lfsu_s.finiteElement().localBasis().order();
566 const int order_n = lfsu_n.finiteElement().localBasis().order();
567 int degree = std::max( order_s, order_n );
570 RF penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
573 for (
const auto& ip : rule)
576 const Dune::FieldVector<DF,dim> n_F_local = ig.unitOuterNormal(ip.position());
579 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(ip.position());
580 Dune::FieldVector<DF,dim> iplocal_n = ig.geometryInOutside().global(ip.position());
584 std::vector<RangeType> phi_s(lfsu_s.size());
585 lfsu_s.finiteElement().localBasis().evaluateFunction(iplocal_s,phi_s);
586 std::vector<RangeType> phi_n(lfsu_n.size());
587 lfsu_n.finiteElement().localBasis().evaluateFunction(iplocal_n,phi_n);
589 const std::vector<RangeType>& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
590 const std::vector<RangeType>& phi_n = cache[order_n].evaluateFunction(iplocal_n,lfsu_n.finiteElement().localBasis());
595 std::vector<JacobianType> gradphi_s(lfsu_s.size());
596 lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
597 std::vector<JacobianType> gradphi_n(lfsu_n.size());
598 lfsu_n.finiteElement().localBasis().evaluateJacobian(iplocal_n,gradphi_n);
600 const std::vector<JacobianType>& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
601 const std::vector<JacobianType>& gradphi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsu_n.finiteElement().localBasis());
605 jac = inside_cell.geometry().jacobianInverseTransposed(iplocal_s);
606 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
607 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
608 jac = outside_cell.geometry().jacobianInverseTransposed(iplocal_n);
609 std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
610 for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
613 typename T::Traits::RangeType b = param.b(inside_cell,iplocal_s);
614 RF normalflux = b*n_F_local;
615 RF omegaup_s, omegaup_n;
628 RF factor = ip.weight() * ig.geometry().integrationElement(ip.position());
629 RF ipfactor = penalty_factor * factor;
632 for (size_type j=0; j<lfsu_s.size(); j++) {
633 RF temp1 = -(An_F_s*tgradphi_s[j])*omega_s*factor;
634 for (size_type i=0; i<lfsu_s.size(); i++) {
635 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,omegaup_s * phi_s[j] * normalflux *factor * phi_s[i]);
636 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,temp1 * phi_s[i]);
637 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * factor * theta * omega_s * (An_F_s*tgradphi_s[i]));
638 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * ipfactor * phi_s[i]);
641 for (size_type j=0; j<lfsu_n.size(); j++) {
642 RF temp1 = -(An_F_n*tgradphi_n[j])*omega_n*factor;
643 for (size_type i=0; i<lfsu_s.size(); i++) {
644 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,omegaup_n * phi_n[j] * normalflux *factor * phi_s[i]);
645 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,temp1 * phi_s[i]);
646 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,-phi_n[j] * factor * theta * omega_s * (An_F_s*tgradphi_s[i]));
647 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,-phi_n[j] * ipfactor * phi_s[i]);
650 for (size_type j=0; j<lfsu_s.size(); j++) {
651 RF temp1 = -(An_F_s*tgradphi_s[j])*omega_s*factor;
652 for (size_type i=0; i<lfsu_n.size(); i++) {
653 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-omegaup_s * phi_s[j] * normalflux *factor * phi_n[i]);
654 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-temp1 * phi_n[i]);
655 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,phi_s[j] * factor * theta * omega_n * (An_F_n*tgradphi_n[i]));
656 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-phi_s[j] * ipfactor * phi_n[i]);
659 for (size_type j=0; j<lfsu_n.size(); j++) {
660 RF temp1 = -(An_F_n*tgradphi_n[j])*omega_n*factor;
661 for (size_type i=0; i<lfsu_n.size(); i++) {
662 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-omegaup_n * phi_n[j] * normalflux *factor * phi_n[i]);
663 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-temp1 * phi_n[i]);
664 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-phi_n[j] * factor * theta * omega_n * (An_F_n*tgradphi_n[i]));
665 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,phi_n[j] * ipfactor * phi_n[i]);
673 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
675 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
679 typedef typename LFSV::Traits::FiniteElementType::
680 Traits::LocalBasisType::Traits::DomainFieldType DF;
681 typedef typename LFSV::Traits::FiniteElementType::
682 Traits::LocalBasisType::Traits::RangeFieldType RF;
683 typedef typename LFSV::Traits::FiniteElementType::
684 Traits::LocalBasisType::Traits::RangeType RangeType;
685 typedef typename LFSU::Traits::FiniteElementType::
686 Traits::LocalBasisType::Traits::JacobianType JacobianType;
687 typedef typename LFSV::Traits::SizeType size_type;
690 const int dim = IG::dimension;
691 const int order = std::max(
692 lfsu_s.finiteElement().localBasis().order(),
693 lfsv_s.finiteElement().localBasis().order()
695 const int intorder = intorderadd+quadrature_factor*order;
698 auto inside_cell = ig.inside();
701 const Dune::FieldVector<DF,dim>&
702 inside_local = Dune::ReferenceElements<DF,dim>::general(inside_cell.type()).position(0,0);
703 typename T::Traits::PermTensorType A_s;
704 A_s = param.A(inside_cell,inside_local);
707 Dune::GeometryType gtface = ig.geometryInInside().type();
708 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
713 element_size(inside_cell.geometry(),h_s,hmax_s);
715 h_F = inside_cell.geometry().volume()/ig.geometry().volume();
718 Dune::FieldMatrix<DF,dim,dim> jac;
721 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
722 Dune::FieldVector<RF,dim> An_F_s;
726 harmonic_average = An_F_s*n_F;
728 harmonic_average = 1.0;
731 const int order_s = lfsu_s.finiteElement().localBasis().order();
732 int degree = order_s;
735 RF penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
738 for (
const auto& ip : rule)
740 BCType bctype = param.bctype(ig.intersection(),ip.position());
746 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(ip.position());
749 const Dune::FieldVector<DF,dim> n_F_local = ig.unitOuterNormal(ip.position());
753 std::vector<RangeType> phi_s(lfsu_s.size());
754 lfsu_s.finiteElement().localBasis().evaluateFunction(iplocal_s,phi_s);
755 std::vector<RangeType> psi_s(lfsv_s.size());
756 lfsv_s.finiteElement().localBasis().evaluateFunction(iplocal_s,psi_s);
758 const std::vector<RangeType>& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
759 const std::vector<RangeType>& psi_s = cache[order_s].evaluateFunction(iplocal_s,lfsv_s.finiteElement().localBasis());
763 RF factor = ip.weight() * ig.geometry().integrationElement(ip.position());
768 RF j = param.j(ig.intersection(),ip.position());
771 for (size_type i=0; i<lfsv_s.size(); i++)
772 r_s.accumulate(lfsv_s,i,j * psi_s[i] * factor);
779 for (size_type i=0; i<lfsu_s.size(); i++)
780 u_s += x_s(lfsu_s,i)*phi_s[i];
783 typename T::Traits::RangeType b = param.b(inside_cell,iplocal_s);
784 RF normalflux = b*n_F_local;
788 if (normalflux<-1
e-30)
789 DUNE_THROW(Dune::Exception,
790 "Outflow boundary condition on inflow! [b("
791 << ig.geometry().global(ip.position()) <<
") = "
795 RF term1 = u_s * normalflux *factor;
796 for (size_type i=0; i<lfsv_s.size(); i++)
797 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
800 RF o = param.o(ig.intersection(),ip.position());
803 for (size_type i=0; i<lfsv_s.size(); i++)
804 r_s.accumulate(lfsv_s,i,o * psi_s[i] * factor);
812 std::vector<JacobianType> gradphi_s(lfsu_s.size());
813 lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
814 std::vector<JacobianType> gradpsi_s(lfsv_s.size());
815 lfsv_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradpsi_s);
817 const std::vector<JacobianType>& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
818 const std::vector<JacobianType>& gradpsi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsv_s.finiteElement().localBasis());
822 jac = inside_cell.geometry().jacobianInverseTransposed(iplocal_s);
823 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
824 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
825 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_s(lfsv_s.size());
826 for (size_type i=0; i<lfsv_s.size(); i++) jac.mv(gradpsi_s[i][0],tgradpsi_s[i]);
829 Dune::FieldVector<RF,dim> gradu_s(0.0);
830 for (size_type i=0; i<lfsu_s.size(); i++)
831 gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
834 RF g = param.g(inside_cell,iplocal_s);
837 RF omegaup_s, omegaup_n;
850 RF term1 = (omegaup_s*u_s + omegaup_n*g) * normalflux *factor;
851 for (size_type i=0; i<lfsv_s.size(); i++)
852 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
855 RF term2 = (An_F_s*gradu_s) * factor;
856 for (size_type i=0; i<lfsv_s.size(); i++)
857 r_s.accumulate(lfsv_s,i,-term2 * psi_s[i]);
860 RF term3 = (u_s-g) * factor;
861 for (size_type i=0; i<lfsv_s.size(); i++)
862 r_s.accumulate(lfsv_s,i,term3 * theta * (An_F_s*tgradpsi_s[i]));
865 RF term4 = penalty_factor * (u_s-g) * factor;
866 for (size_type i=0; i<lfsv_s.size(); i++)
867 r_s.accumulate(lfsv_s,i,term4 * psi_s[i]);
871 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
873 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
877 typedef typename LFSV::Traits::FiniteElementType::
878 Traits::LocalBasisType::Traits::DomainFieldType DF;
879 typedef typename LFSV::Traits::FiniteElementType::
880 Traits::LocalBasisType::Traits::RangeFieldType RF;
881 typedef typename LFSV::Traits::FiniteElementType::
882 Traits::LocalBasisType::Traits::RangeType RangeType;
883 typedef typename LFSU::Traits::FiniteElementType::
884 Traits::LocalBasisType::Traits::JacobianType JacobianType;
885 typedef typename LFSV::Traits::SizeType size_type;
888 const int dim = IG::dimension;
889 const int order = std::max(
890 lfsu_s.finiteElement().localBasis().order(),
891 lfsv_s.finiteElement().localBasis().order()
893 const int intorder = intorderadd+quadrature_factor*order;
896 auto inside_cell = ig.inside();
899 const Dune::FieldVector<DF,dim>&
900 inside_local = Dune::ReferenceElements<DF,dim>::general(inside_cell.type()).position(0,0);
901 typename T::Traits::PermTensorType A_s;
902 A_s = param.A(inside_cell,inside_local);
905 Dune::GeometryType gtface = ig.geometryInInside().type();
906 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
911 element_size(inside_cell.geometry(),h_s,hmax_s);
913 h_F = inside_cell.geometry().volume()/ig.geometry().volume();
916 Dune::FieldMatrix<DF,dim,dim> jac;
919 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
920 Dune::FieldVector<RF,dim> An_F_s;
924 harmonic_average = An_F_s*n_F;
926 harmonic_average = 1.0;
929 const int order_s = lfsu_s.finiteElement().localBasis().order();
930 int degree = order_s;
933 RF penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
939 for (
const auto& ip : rule)
941 BCType bctype = param.bctype(ig.intersection(),ip.position());
948 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(ip.position());
951 const Dune::FieldVector<DF,dim> n_F_local = ig.unitOuterNormal(ip.position());
955 std::vector<RangeType> phi_s(lfsu_s.size());
956 lfsu_s.finiteElement().localBasis().evaluateFunction(iplocal_s,phi_s);
958 const std::vector<RangeType>& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
962 RF factor = ip.weight() * ig.geometry().integrationElement(ip.position());
965 typename T::Traits::RangeType b = param.b(inside_cell,iplocal_s);
966 RF normalflux = b*n_F_local;
970 if (normalflux<-1
e-30)
971 DUNE_THROW(Dune::Exception,
972 "Outflow boundary condition on inflow! [b("
973 << ig.geometry().global(ip.position()) <<
") = "
974 << b <<
")" << n_F_local <<
" " << normalflux);
977 for (size_type j=0; j<lfsu_s.size(); j++)
978 for (size_type i=0; i<lfsu_s.size(); i++)
979 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * normalflux * factor * phi_s[i]);
986 std::vector<JacobianType> gradphi_s(lfsu_s.size());
987 lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
989 const std::vector<JacobianType>& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
993 jac = inside_cell.geometry().jacobianInverseTransposed(iplocal_s);
994 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
995 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
998 RF omegaup_s, omegaup_n;
1011 for (size_type j=0; j<lfsu_s.size(); j++)
1012 for (size_type i=0; i<lfsu_s.size(); i++)
1013 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,omegaup_s * phi_s[j] * normalflux * factor * phi_s[i]);
1016 for (size_type j=0; j<lfsu_s.size(); j++)
1017 for (size_type i=0; i<lfsu_s.size(); i++)
1018 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,-(An_F_s*tgradphi_s[j]) * factor * phi_s[i]);
1021 for (size_type j=0; j<lfsu_s.size(); j++)
1022 for (size_type i=0; i<lfsu_s.size(); i++)
1023 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * factor * theta * (An_F_s*tgradphi_s[i]));
1026 for (size_type j=0; j<lfsu_s.size(); j++)
1027 for (size_type i=0; i<lfsu_s.size(); i++)
1028 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,penalty_factor * phi_s[j] * phi_s[i] * factor);
1033 template<
typename EG,
typename LFSV,
typename R>
1037 typedef typename LFSV::Traits::FiniteElementType::
1038 Traits::LocalBasisType::Traits::DomainFieldType DF;
1039 typedef typename LFSV::Traits::FiniteElementType::
1040 Traits::LocalBasisType::Traits::RangeFieldType RF;
1041 typedef typename LFSV::Traits::FiniteElementType::
1042 Traits::LocalBasisType::Traits::RangeType RangeType;
1043 typedef typename LFSV::Traits::SizeType size_type;
1046 const int dim = EG::Entity::dimension;
1047 const int order = lfsv.finiteElement().localBasis().order();
1048 const int intorder = intorderadd + 2 * order;
1051 Dune::GeometryType gt = eg.geometry().type();
1052 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
1055 for (
const auto& ip : rule)
1059 std::vector<RangeType> phi(lfsv.size());
1060 lfsv.finiteElement().localBasis().evaluateFunction(ip.position(),phi);
1062 const std::vector<RangeType>& phi = cache[order].evaluateFunction(ip.position(),lfsv.finiteElement().localBasis());
1067 f = param.f(eg.entity(),ip.position());
1070 RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
1071 for (size_type i=0; i<lfsv.size(); i++)
1072 r.accumulate(lfsv,i,-f*phi[i]*factor);
1089 int quadrature_factor;
1091 typedef typename FiniteElementMap::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
1104 std::vector<Cache> cache;
1107 void element_size (
const GEO& geo,
typename GEO::ctype& hmin,
typename GEO::ctype hmax)
const
1109 typedef typename GEO::ctype DF;
1112 const int dim = GEO::coorddimension;
1115 Dune::FieldVector<DF,dim> x = geo.corner(0);
1117 hmin = hmax = x.two_norm();
1122 Dune::GeometryType gt = geo.type();
1123 for (
int i=0; i<Dune::ReferenceElements<DF,dim>::general(gt).size(dim-1); i++)
1125 Dune::FieldVector<DF,dim> x = geo.corner(Dune::ReferenceElements<DF,dim>::general(gt).subEntity(i,dim-1,0,dim));
1126 x -= geo.corner(Dune::ReferenceElements<DF,dim>::general(gt).subEntity(i,dim-1,1,dim));
1127 hmin = std::min(hmin,x.two_norm());
1128 hmax = std::max(hmax,x.two_norm());
Definition: convectiondiffusiondg.hh:77
Definition: convectiondiffusiondg.hh:74
Definition: convectiondiffusiondg.hh:34
void alpha_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, R &r_s) const
Definition: convectiondiffusiondg.hh:674
void setTime(R t_)
set time for subsequent evaluation
Definition: idefault.hh:104
ConvectionDiffusionDG(T ¶m_, ConvectionDiffusionDGMethod::Type method_=ConvectionDiffusionDGMethod::NIPG, ConvectionDiffusionDGWeights::Type weights_=ConvectionDiffusionDGWeights::weightsOff, Real alpha_=0.0, int intorderadd_=0)
constructor: pass parameter object
Definition: convectiondiffusiondg.hh:83
Definition: convectiondiffusiondg.hh:39
const E & e
Definition: interpolate.hh:172
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: convectiondiffusiondg.hh:1034
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusiondg.hh:102
Definition: convectiondiffusionparameter.hh:67
Type
Definition: convectiondiffusiondg.hh:34
Definition: convectiondiffusiondg.hh:73
Definition: convectiondiffusiondg.hh:80
const IG & ig
Definition: constraints.hh:147
Type
Definition: convectiondiffusionparameter.hh:67
Default flags for all local operators.
Definition: flags.hh:18
Definition: convectiondiffusionparameter.hh:67
sparsity pattern generator
Definition: pattern.hh:13
Definition: convectiondiffusiondg.hh:57
Definition: convectiondiffusionparameter.hh:67
Definition: convectiondiffusiondg.hh:78
Definition: convectiondiffusionparameter.hh:67
Definition: convectiondiffusiondg.hh:37
void setTime(double t)
set time in parameter class
Definition: convectiondiffusiondg.hh:1077
void alpha_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, R &r_s, R &r_n) const
Definition: convectiondiffusiondg.hh:275
Implement jacobian_apply_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:445
Definition: adaptivity.hh:27
Definition: convectiondiffusiondg.hh:39
Definition: convectiondiffusiondg.hh:79
sparsity pattern generator
Definition: pattern.hh:29
Type
Definition: convectiondiffusiondg.hh:39
Definition: convectiondiffusiondg.hh:34
void jacobian_skeleton(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, const LFSU &lfsu_n, const X &x_n, const LFSV &lfsv_n, M &mat_ss, M &mat_sn, M &mat_ns, M &mat_nn) const
Definition: convectiondiffusiondg.hh:484
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:15
Definition: convectiondiffusiondg.hh:32
const EG & eg
Definition: constraints.hh:280
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusiondg.hh:197
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538
void jacobian_boundary(const IG &ig, const LFSU &lfsu_s, const X &x_s, const LFSV &lfsv_s, M &mat_ss) const
Definition: convectiondiffusiondg.hh:872