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/common/static_assert.hh>
8 #include<dune/geometry/referenceelements.hh>
57 template<
typename T,
typename FiniteElementMap>
67 enum { dim = T::Traits::GridViewType::dimension };
69 typedef typename T::Traits::RangeFieldType Real;
93 param(param_), method(method_), weights(weights_),
94 alpha(alpha_), intorderadd(intorderadd_), quadrature_factor(2),
102 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
103 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
106 typedef typename LFSU::Traits::FiniteElementType::
107 Traits::LocalBasisType::Traits::DomainFieldType DF;
108 typedef typename LFSU::Traits::FiniteElementType::
109 Traits::LocalBasisType::Traits::RangeFieldType RF;
110 typedef typename LFSU::Traits::FiniteElementType::
111 Traits::LocalBasisType::Traits::JacobianType JacobianType;
112 typedef typename LFSU::Traits::FiniteElementType::
113 Traits::LocalBasisType::Traits::RangeType RangeType;
114 typedef typename LFSU::Traits::SizeType size_type;
117 const int dim = EG::Geometry::dimension;
118 const int order = std::max(lfsu.finiteElement().localBasis().order(),
119 lfsv.finiteElement().localBasis().order());
120 const int intorder = intorderadd + quadrature_factor * order;
123 Dune::GeometryType gt = eg.geometry().type();
124 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
127 typename T::Traits::PermTensorType A;
128 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
129 A = param.A(eg.entity(),localcenter);
132 typename EG::Geometry::JacobianInverseTransposed jac;
135 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
139 std::vector<RangeType> phi(lfsu.size());
140 lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
141 std::vector<RangeType> psi(lfsv.size());
142 lfsv.finiteElement().localBasis().evaluateFunction(it->position(),psi);
144 const std::vector<RangeType>& phi = cache[order].evaluateFunction(it->position(),lfsu.finiteElement().localBasis());
145 const std::vector<RangeType>& psi = cache[order].evaluateFunction(it->position(),lfsv.finiteElement().localBasis());
150 for (size_type i=0; i<lfsu.size(); i++)
151 u += x(lfsu,i)*phi[i];
155 std::vector<JacobianType> js(lfsu.size());
156 lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
157 std::vector<JacobianType> js_v(lfsv.size());
158 lfsv.finiteElement().localBasis().evaluateJacobian(it->position(),js_v);
160 const std::vector<JacobianType>& js = cache[order].evaluateJacobian(it->position(),lfsu.finiteElement().localBasis());
161 const std::vector<JacobianType>& js_v = cache[order].evaluateJacobian(it->position(),lfsv.finiteElement().localBasis());
165 jac = eg.geometry().jacobianInverseTransposed(it->position());
166 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
167 for (size_type i=0; i<lfsu.size(); i++)
168 jac.mv(js[i][0],gradphi[i]);
170 std::vector<Dune::FieldVector<RF,dim> > gradpsi(lfsv.size());
171 for (size_type i=0; i<lfsv.size(); i++)
172 jac.mv(js_v[i][0],gradpsi[i]);
175 Dune::FieldVector<RF,dim> gradu(0.0);
176 for (size_type i=0; i<lfsu.size(); i++)
177 gradu.axpy(x(lfsu,i),gradphi[i]);
180 Dune::FieldVector<RF,dim> Agradu(0.0);
184 typename T::Traits::RangeType b = param.b(eg.entity(),it->position());
187 typename T::Traits::RangeFieldType c = param.c(eg.entity(),it->position());
190 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
191 for (size_type i=0; i<lfsv.size(); i++)
192 r.accumulate(lfsv,i,( Agradu*gradpsi[i] - u*(b*gradpsi[i]) + c*u*psi[i] )*factor);
197 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
202 typedef typename LFSU::Traits::FiniteElementType::
203 Traits::LocalBasisType::Traits::DomainFieldType DF;
204 typedef typename LFSU::Traits::FiniteElementType::
205 Traits::LocalBasisType::Traits::RangeFieldType RF;
206 typedef typename LFSU::Traits::FiniteElementType::
207 Traits::LocalBasisType::Traits::JacobianType JacobianType;
208 typedef typename LFSU::Traits::FiniteElementType::
209 Traits::LocalBasisType::Traits::RangeType RangeType;
210 typedef typename LFSU::Traits::SizeType size_type;
213 const int dim = EG::Geometry::dimension;
214 const int order = std::max(lfsu.finiteElement().localBasis().order(),
215 lfsv.finiteElement().localBasis().order());
216 const int intorder = intorderadd + quadrature_factor * order;
219 Dune::GeometryType gt = eg.geometry().type();
220 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
223 typename T::Traits::PermTensorType A;
224 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
225 A = param.A(eg.entity(),localcenter);
228 typename EG::Geometry::JacobianInverseTransposed jac;
231 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
235 std::vector<RangeType> phi(lfsu.size());
236 lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
238 const std::vector<RangeType>& phi = cache[order].evaluateFunction(it->position(),lfsu.finiteElement().localBasis());
243 std::vector<JacobianType> js(lfsu.size());
244 lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
246 const std::vector<JacobianType>& js = cache[order].evaluateJacobian(it->position(),lfsu.finiteElement().localBasis());
250 jac = eg.geometry().jacobianInverseTransposed(it->position());
251 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
252 std::vector<Dune::FieldVector<RF,dim> > Agradphi(lfsu.size());
253 for (size_type i=0; i<lfsu.size(); i++)
255 jac.mv(js[i][0],gradphi[i]);
256 A.mv(gradphi[i],Agradphi[i]);
260 typename T::Traits::RangeType b = param.b(eg.entity(),it->position());
263 typename T::Traits::RangeFieldType c = param.c(eg.entity(),it->position());
266 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
267 for (size_type j=0; j<lfsu.size(); j++)
268 for (size_type i=0; i<lfsu.size(); i++)
269 mat.accumulate(lfsu,i,lfsu,j,( Agradphi[j]*gradphi[i] - phi[j]*(b*gradphi[i]) + c*phi[j]*phi[i] )*factor);
275 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
277 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
278 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
279 R& r_s, R& r_n)
const
282 typedef typename LFSV::Traits::FiniteElementType::
283 Traits::LocalBasisType::Traits::DomainFieldType DF;
284 typedef typename LFSV::Traits::FiniteElementType::
285 Traits::LocalBasisType::Traits::RangeFieldType RF;
286 typedef typename LFSV::Traits::FiniteElementType::
287 Traits::LocalBasisType::Traits::RangeType RangeType;
288 typedef typename LFSU::Traits::FiniteElementType::
289 Traits::LocalBasisType::Traits::JacobianType JacobianType;
290 typedef typename LFSV::Traits::SizeType size_type;
293 const int dim = IG::dimension;
294 const int order = std::max(
295 std::max(lfsu_s.finiteElement().localBasis().order(),
296 lfsu_n.finiteElement().localBasis().order()),
297 std::max(lfsv_s.finiteElement().localBasis().order(),
298 lfsv_n.finiteElement().localBasis().order())
300 const int intorder = intorderadd+quadrature_factor*order;
303 const Dune::FieldVector<DF,dim>&
304 inside_local = Dune::ReferenceElements<DF,dim>::general(ig.inside()->type()).position(0,0);
305 const Dune::FieldVector<DF,dim>&
306 outside_local = Dune::ReferenceElements<DF,dim>::general(ig.outside()->type()).position(0,0);
307 typename T::Traits::PermTensorType A_s, A_n;
308 A_s = param.A(*(ig.inside()),inside_local);
309 A_n = param.A(*(ig.outside()),outside_local);
315 element_size(ig.inside()->geometry(),h_s,hmax_s);
316 element_size(ig.outside()->geometry(),h_n,hmax_n);
317 RF h_F = std::min(h_s,h_n);
318 h_F = std::min(ig.inside()->geometry().volume(),ig.outside()->geometry().volume())/ig.geometry().volume();
321 Dune::GeometryType gtface = ig.geometryInInside().type();
322 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
325 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
328 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
329 Dune::FieldVector<RF,dim> An_F_s;
331 Dune::FieldVector<RF,dim> An_F_n;
337 RF harmonic_average(0.0);
340 RF delta_s = (An_F_s*n_F);
341 RF delta_n = (An_F_n*n_F);
342 omega_s = delta_n/(delta_s+delta_n+1
e-20);
343 omega_n = delta_s/(delta_s+delta_n+1
e-20);
344 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1
e-20);
348 omega_s = omega_n = 0.5;
349 harmonic_average = 1.0;
353 const int order_s = lfsu_s.finiteElement().localBasis().order();
354 const int order_n = lfsu_n.finiteElement().localBasis().order();
355 int degree = std::max( order_s, order_n );
358 RF penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
361 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
364 const Dune::FieldVector<DF,dim> n_F_local = ig.unitOuterNormal(it->position());
367 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(it->position());
368 Dune::FieldVector<DF,dim> iplocal_n = ig.geometryInOutside().global(it->position());
372 std::vector<RangeType> phi_s(lfsu_s.size());
373 lfsu_s.finiteElement().localBasis().evaluateFunction(iplocal_s,phi_s);
374 std::vector<RangeType> phi_n(lfsu_n.size());
375 lfsu_n.finiteElement().localBasis().evaluateFunction(iplocal_n,phi_n);
376 std::vector<RangeType> psi_s(lfsv_s.size());
377 lfsv_s.finiteElement().localBasis().evaluateFunction(iplocal_s,psi_s);
378 std::vector<RangeType> psi_n(lfsv_n.size());
379 lfsv_n.finiteElement().localBasis().evaluateFunction(iplocal_n,psi_n);
381 const std::vector<RangeType>& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
382 const std::vector<RangeType>& phi_n = cache[order_n].evaluateFunction(iplocal_n,lfsu_n.finiteElement().localBasis());
383 const std::vector<RangeType>& psi_s = cache[order_s].evaluateFunction(iplocal_s,lfsv_s.finiteElement().localBasis());
384 const std::vector<RangeType>& psi_n = cache[order_n].evaluateFunction(iplocal_n,lfsv_n.finiteElement().localBasis());
389 for (size_type i=0; i<lfsu_s.size(); i++)
390 u_s += x_s(lfsu_s,i)*phi_s[i];
392 for (size_type i=0; i<lfsu_n.size(); i++)
393 u_n += x_n(lfsu_n,i)*phi_n[i];
397 std::vector<JacobianType> gradphi_s(lfsu_s.size());
398 lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
399 std::vector<JacobianType> gradphi_n(lfsu_n.size());
400 lfsu_n.finiteElement().localBasis().evaluateJacobian(iplocal_n,gradphi_n);
401 std::vector<JacobianType> gradpsi_s(lfsv_s.size());
402 lfsv_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradpsi_s);
403 std::vector<JacobianType> gradpsi_n(lfsv_n.size());
404 lfsv_n.finiteElement().localBasis().evaluateJacobian(iplocal_n,gradpsi_n);
406 const std::vector<JacobianType>& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
407 const std::vector<JacobianType>& gradphi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsu_n.finiteElement().localBasis());
408 const std::vector<JacobianType>& gradpsi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsv_s.finiteElement().localBasis());
409 const std::vector<JacobianType>& gradpsi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsv_n.finiteElement().localBasis());
413 jac = ig.inside()->geometry().jacobianInverseTransposed(iplocal_s);
414 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
415 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
416 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_s(lfsv_s.size());
417 for (size_type i=0; i<lfsv_s.size(); i++) jac.mv(gradpsi_s[i][0],tgradpsi_s[i]);
418 jac = ig.outside()->geometry().jacobianInverseTransposed(iplocal_n);
419 std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
420 for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
421 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_n(lfsv_n.size());
422 for (size_type i=0; i<lfsv_n.size(); i++) jac.mv(gradpsi_n[i][0],tgradpsi_n[i]);
425 Dune::FieldVector<RF,dim> gradu_s(0.0);
426 for (size_type i=0; i<lfsu_s.size(); i++)
427 gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
428 Dune::FieldVector<RF,dim> gradu_n(0.0);
429 for (size_type i=0; i<lfsu_n.size(); i++)
430 gradu_n.axpy(x_n(lfsu_n,i),tgradphi_n[i]);
433 typename T::Traits::RangeType b = param.b(*(ig.inside()),iplocal_s);
434 RF normalflux = b*n_F_local;
435 RF omegaup_s, omegaup_n;
448 RF factor = it->weight() * ig.geometry().integrationElement(it->position());
451 RF term1 = (omegaup_s*u_s + omegaup_n*u_n) * normalflux *factor;
452 for (size_type i=0; i<lfsv_s.size(); i++)
453 r_s.accumulate(lfsu_s,i,term1 * psi_s[i]);
454 for (size_type i=0; i<lfsv_n.size(); i++)
455 r_n.accumulate(lfsu_n,i,-term1 * psi_n[i]);
458 RF term2 = -(omega_s*(An_F_s*gradu_s) + omega_n*(An_F_n*gradu_n)) * factor;
459 for (size_type i=0; i<lfsv_s.size(); i++)
460 r_s.accumulate(lfsv_s,i,term2 * psi_s[i]);
461 for (size_type i=0; i<lfsv_n.size(); i++)
462 r_n.accumulate(lfsv_n,i,-term2 * psi_n[i]);
465 RF term3 = (u_s-u_n) * factor;
466 for (size_type i=0; i<lfsv_s.size(); i++)
467 r_s.accumulate(lfsv_s,i,term3 * theta * omega_s * (An_F_s*tgradpsi_s[i]));
468 for (size_type i=0; i<lfsv_n.size(); i++)
469 r_n.accumulate(lfsv_n,i,term3 * theta * omega_n * (An_F_n*tgradpsi_n[i]));
472 RF term4 = penalty_factor * (u_s-u_n) * factor;
473 for (size_type i=0; i<lfsv_s.size(); i++)
474 r_s.accumulate(lfsv_s,i,term4 * psi_s[i]);
475 for (size_type i=0; i<lfsv_n.size(); i++)
476 r_n.accumulate(lfsv_n,i,-term4 * psi_n[i]);
480 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
482 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
483 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
484 M& mat_ss, M& mat_sn,
485 M& mat_ns, M& mat_nn)
const
488 typedef typename LFSV::Traits::FiniteElementType::
489 Traits::LocalBasisType::Traits::DomainFieldType DF;
490 typedef typename LFSV::Traits::FiniteElementType::
491 Traits::LocalBasisType::Traits::RangeFieldType RF;
492 typedef typename LFSV::Traits::FiniteElementType::
493 Traits::LocalBasisType::Traits::RangeType RangeType;
494 typedef typename LFSU::Traits::FiniteElementType::
495 Traits::LocalBasisType::Traits::JacobianType JacobianType;
496 typedef typename LFSV::Traits::SizeType size_type;
499 const int dim = IG::dimension;
500 const int order = std::max(
501 std::max(lfsu_s.finiteElement().localBasis().order(),
502 lfsu_n.finiteElement().localBasis().order()),
503 std::max(lfsv_s.finiteElement().localBasis().order(),
504 lfsv_n.finiteElement().localBasis().order())
506 const int intorder = intorderadd+quadrature_factor*order;
509 const Dune::FieldVector<DF,dim>&
510 inside_local = Dune::ReferenceElements<DF,dim>::general(ig.inside()->type()).position(0,0);
511 const Dune::FieldVector<DF,dim>&
512 outside_local = Dune::ReferenceElements<DF,dim>::general(ig.outside()->type()).position(0,0);
513 typename T::Traits::PermTensorType A_s, A_n;
514 A_s = param.A(*(ig.inside()),inside_local);
515 A_n = param.A(*(ig.outside()),outside_local);
519 DF hmax_s = 0., hmax_n = 0.;
520 element_size(ig.inside()->geometry(),h_s,hmax_s);
521 element_size(ig.outside()->geometry(),h_n,hmax_n);
522 RF h_F = std::min(h_s,h_n);
523 h_F = std::min(ig.inside()->geometry().volume(),ig.outside()->geometry().volume())/ig.geometry().volume();
526 Dune::GeometryType gtface = ig.geometryInInside().type();
527 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
530 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
533 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
534 Dune::FieldVector<RF,dim> An_F_s;
536 Dune::FieldVector<RF,dim> An_F_n;
542 RF harmonic_average(0.0);
545 RF delta_s = (An_F_s*n_F);
546 RF delta_n = (An_F_n*n_F);
547 omega_s = delta_n/(delta_s+delta_n+1
e-20);
548 omega_n = delta_s/(delta_s+delta_n+1
e-20);
549 harmonic_average = 2.0*delta_s*delta_n/(delta_s+delta_n+1
e-20);
553 omega_s = omega_n = 0.5;
554 harmonic_average = 1.0;
558 const int order_s = lfsu_s.finiteElement().localBasis().order();
559 const int order_n = lfsu_n.finiteElement().localBasis().order();
560 int degree = std::max( order_s, order_n );
563 RF penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
566 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
569 const Dune::FieldVector<DF,dim> n_F_local = ig.unitOuterNormal(it->position());
572 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(it->position());
573 Dune::FieldVector<DF,dim> iplocal_n = ig.geometryInOutside().global(it->position());
577 std::vector<RangeType> phi_s(lfsu_s.size());
578 lfsu_s.finiteElement().localBasis().evaluateFunction(iplocal_s,phi_s);
579 std::vector<RangeType> phi_n(lfsu_n.size());
580 lfsu_n.finiteElement().localBasis().evaluateFunction(iplocal_n,phi_n);
582 const std::vector<RangeType>& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
583 const std::vector<RangeType>& phi_n = cache[order_n].evaluateFunction(iplocal_n,lfsu_n.finiteElement().localBasis());
588 std::vector<JacobianType> gradphi_s(lfsu_s.size());
589 lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
590 std::vector<JacobianType> gradphi_n(lfsu_n.size());
591 lfsu_n.finiteElement().localBasis().evaluateJacobian(iplocal_n,gradphi_n);
593 const std::vector<JacobianType>& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
594 const std::vector<JacobianType>& gradphi_n = cache[order_n].evaluateJacobian(iplocal_n,lfsu_n.finiteElement().localBasis());
598 jac = ig.inside()->geometry().jacobianInverseTransposed(iplocal_s);
599 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
600 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
601 jac = ig.outside()->geometry().jacobianInverseTransposed(iplocal_n);
602 std::vector<Dune::FieldVector<RF,dim> > tgradphi_n(lfsu_n.size());
603 for (size_type i=0; i<lfsu_n.size(); i++) jac.mv(gradphi_n[i][0],tgradphi_n[i]);
606 typename T::Traits::RangeType b = param.b(*(ig.inside()),iplocal_s);
607 RF normalflux = b*n_F_local;
608 RF omegaup_s, omegaup_n;
621 RF factor = it->weight() * ig.geometry().integrationElement(it->position());
622 RF ipfactor = penalty_factor * factor;
625 for (size_type j=0; j<lfsu_s.size(); j++) {
626 RF temp1 = -(An_F_s*tgradphi_s[j])*omega_s*factor;
627 for (size_type i=0; i<lfsu_s.size(); i++) {
628 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,omegaup_s * phi_s[j] * normalflux *factor * phi_s[i]);
629 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,temp1 * phi_s[i]);
630 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * factor * theta * omega_s * (An_F_s*tgradphi_s[i]));
631 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * ipfactor * phi_s[i]);
634 for (size_type j=0; j<lfsu_n.size(); j++) {
635 RF temp1 = -(An_F_n*tgradphi_n[j])*omega_n*factor;
636 for (size_type i=0; i<lfsu_s.size(); i++) {
637 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,omegaup_n * phi_n[j] * normalflux *factor * phi_s[i]);
638 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,temp1 * phi_s[i]);
639 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,-phi_n[j] * factor * theta * omega_s * (An_F_s*tgradphi_s[i]));
640 mat_sn.accumulate(lfsu_s,i,lfsu_n,j,-phi_n[j] * ipfactor * phi_s[i]);
643 for (size_type j=0; j<lfsu_s.size(); j++) {
644 RF temp1 = -(An_F_s*tgradphi_s[j])*omega_s*factor;
645 for (size_type i=0; i<lfsu_n.size(); i++) {
646 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-omegaup_s * phi_s[j] * normalflux *factor * phi_n[i]);
647 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-temp1 * phi_n[i]);
648 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,phi_s[j] * factor * theta * omega_n * (An_F_n*tgradphi_n[i]));
649 mat_ns.accumulate(lfsu_n,i,lfsu_s,j,-phi_s[j] * ipfactor * phi_n[i]);
652 for (size_type j=0; j<lfsu_n.size(); j++) {
653 RF temp1 = -(An_F_n*tgradphi_n[j])*omega_n*factor;
654 for (size_type i=0; i<lfsu_n.size(); i++) {
655 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-omegaup_n * phi_n[j] * normalflux *factor * phi_n[i]);
656 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-temp1 * phi_n[i]);
657 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,-phi_n[j] * factor * theta * omega_n * (An_F_n*tgradphi_n[i]));
658 mat_nn.accumulate(lfsu_n,i,lfsu_n,j,phi_n[j] * ipfactor * phi_n[i]);
666 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
668 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
672 typedef typename LFSV::Traits::FiniteElementType::
673 Traits::LocalBasisType::Traits::DomainFieldType DF;
674 typedef typename LFSV::Traits::FiniteElementType::
675 Traits::LocalBasisType::Traits::RangeFieldType RF;
676 typedef typename LFSV::Traits::FiniteElementType::
677 Traits::LocalBasisType::Traits::RangeType RangeType;
678 typedef typename LFSU::Traits::FiniteElementType::
679 Traits::LocalBasisType::Traits::JacobianType JacobianType;
680 typedef typename LFSV::Traits::SizeType size_type;
683 const int dim = IG::dimension;
684 const int order = std::max(
685 lfsu_s.finiteElement().localBasis().order(),
686 lfsv_s.finiteElement().localBasis().order()
688 const int intorder = intorderadd+quadrature_factor*order;
691 const Dune::FieldVector<DF,dim>&
692 inside_local = Dune::ReferenceElements<DF,dim>::general(ig.inside()->type()).position(0,0);
693 typename T::Traits::PermTensorType A_s;
694 A_s = param.A(*(ig.inside()),inside_local);
697 Dune::GeometryType gtface = ig.geometryInInside().type();
698 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
703 element_size(ig.inside()->geometry(),h_s,hmax_s);
705 h_F = ig.inside()->geometry().volume()/ig.geometry().volume();
708 Dune::FieldMatrix<DF,dim,dim> jac;
711 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
712 Dune::FieldVector<RF,dim> An_F_s;
716 harmonic_average = An_F_s*n_F;
718 harmonic_average = 1.0;
721 const int order_s = lfsu_s.finiteElement().localBasis().order();
722 int degree = order_s;
725 RF penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
728 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
730 BCType bctype = param.bctype(ig.intersection(),it->position());
733 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(it->position());
736 const Dune::FieldVector<DF,dim> n_F_local = ig.unitOuterNormal(it->position());
740 std::vector<RangeType> phi_s(lfsu_s.size());
741 lfsu_s.finiteElement().localBasis().evaluateFunction(iplocal_s,phi_s);
742 std::vector<RangeType> psi_s(lfsv_s.size());
743 lfsv_s.finiteElement().localBasis().evaluateFunction(iplocal_s,psi_s);
745 const std::vector<RangeType>& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
746 const std::vector<RangeType>& psi_s = cache[order_s].evaluateFunction(iplocal_s,lfsv_s.finiteElement().localBasis());
750 RF factor = it->weight() * ig.geometry().integrationElement(it->position());
755 RF j = param.j(ig.intersection(),it->position());
758 for (size_type i=0; i<lfsv_s.size(); i++)
759 r_s.accumulate(lfsv_s,i,j * psi_s[i] * factor);
766 for (size_type i=0; i<lfsu_s.size(); i++)
767 u_s += x_s(lfsu_s,i)*phi_s[i];
770 typename T::Traits::RangeType b = param.b(*(ig.inside()),iplocal_s);
771 RF normalflux = b*n_F_local;
775 if (normalflux<-1
e-30)
776 DUNE_THROW(Dune::Exception,
777 "Outflow boundary condition on inflow! [b("
778 << ig.geometry().global(it->position()) <<
") = "
782 RF term1 = u_s * normalflux *factor;
783 for (size_type i=0; i<lfsv_s.size(); i++)
784 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
787 RF o = param.o(ig.intersection(),it->position());
790 for (size_type i=0; i<lfsv_s.size(); i++)
791 r_s.accumulate(lfsv_s,i,o * psi_s[i] * factor);
799 std::vector<JacobianType> gradphi_s(lfsu_s.size());
800 lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
801 std::vector<JacobianType> gradpsi_s(lfsv_s.size());
802 lfsv_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradpsi_s);
804 const std::vector<JacobianType>& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
805 const std::vector<JacobianType>& gradpsi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsv_s.finiteElement().localBasis());
809 jac = ig.inside()->geometry().jacobianInverseTransposed(iplocal_s);
810 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
811 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
812 std::vector<Dune::FieldVector<RF,dim> > tgradpsi_s(lfsv_s.size());
813 for (size_type i=0; i<lfsv_s.size(); i++) jac.mv(gradpsi_s[i][0],tgradpsi_s[i]);
816 Dune::FieldVector<RF,dim> gradu_s(0.0);
817 for (size_type i=0; i<lfsu_s.size(); i++)
818 gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
821 RF g = param.g(*(ig.inside()),iplocal_s);
824 RF omegaup_s, omegaup_n;
837 RF term1 = (omegaup_s*u_s + omegaup_n*g) * normalflux *factor;
838 for (size_type i=0; i<lfsv_s.size(); i++)
839 r_s.accumulate(lfsv_s,i,term1 * psi_s[i]);
842 RF term2 = (An_F_s*gradu_s) * factor;
843 for (size_type i=0; i<lfsv_s.size(); i++)
844 r_s.accumulate(lfsv_s,i,-term2 * psi_s[i]);
847 RF term3 = (u_s-g) * factor;
848 for (size_type i=0; i<lfsv_s.size(); i++)
849 r_s.accumulate(lfsv_s,i,term3 * theta * (An_F_s*tgradpsi_s[i]));
852 RF term4 = penalty_factor * (u_s-g) * factor;
853 for (size_type i=0; i<lfsv_s.size(); i++)
854 r_s.accumulate(lfsv_s,i,term4 * psi_s[i]);
858 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
860 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
864 typedef typename LFSV::Traits::FiniteElementType::
865 Traits::LocalBasisType::Traits::DomainFieldType DF;
866 typedef typename LFSV::Traits::FiniteElementType::
867 Traits::LocalBasisType::Traits::RangeFieldType RF;
868 typedef typename LFSV::Traits::FiniteElementType::
869 Traits::LocalBasisType::Traits::RangeType RangeType;
870 typedef typename LFSU::Traits::FiniteElementType::
871 Traits::LocalBasisType::Traits::JacobianType JacobianType;
872 typedef typename LFSV::Traits::SizeType size_type;
875 const int dim = IG::dimension;
876 const int order = std::max(
877 lfsu_s.finiteElement().localBasis().order(),
878 lfsv_s.finiteElement().localBasis().order()
880 const int intorder = intorderadd+quadrature_factor*order;
883 const Dune::FieldVector<DF,dim>&
884 inside_local = Dune::ReferenceElements<DF,dim>::general(ig.inside()->type()).position(0,0);
885 typename T::Traits::PermTensorType A_s;
886 A_s = param.A(*(ig.inside()),inside_local);
889 Dune::GeometryType gtface = ig.geometryInInside().type();
890 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
895 element_size(ig.inside()->geometry(),h_s,hmax_s);
897 h_F = ig.inside()->geometry().volume()/ig.geometry().volume();
900 Dune::FieldMatrix<DF,dim,dim> jac;
903 const Dune::FieldVector<DF,dim> n_F = ig.centerUnitOuterNormal();
904 Dune::FieldVector<RF,dim> An_F_s;
908 harmonic_average = An_F_s*n_F;
910 harmonic_average = 1.0;
913 const int order_s = lfsu_s.finiteElement().localBasis().order();
914 int degree = order_s;
917 RF penalty_factor = (alpha/h_F) * harmonic_average * degree*(degree+dim-1);
923 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
925 BCType bctype = param.bctype(ig.intersection(),it->position());
929 Dune::FieldVector<DF,dim> iplocal_s = ig.geometryInInside().global(it->position());
932 const Dune::FieldVector<DF,dim> n_F_local = ig.unitOuterNormal(it->position());
936 std::vector<RangeType> phi_s(lfsu_s.size());
937 lfsu_s.finiteElement().localBasis().evaluateFunction(iplocal_s,phi_s);
939 const std::vector<RangeType>& phi_s = cache[order_s].evaluateFunction(iplocal_s,lfsu_s.finiteElement().localBasis());
943 RF factor = it->weight() * ig.geometry().integrationElement(it->position());
946 typename T::Traits::RangeType b = param.b(*(ig.inside()),iplocal_s);
947 RF normalflux = b*n_F_local;
951 if (normalflux<-1
e-30)
952 DUNE_THROW(Dune::Exception,
953 "Outflow boundary condition on inflow! [b("
954 << ig.geometry().global(it->position()) <<
") = "
955 << b <<
")" << n_F_local <<
" " << normalflux);
958 for (size_type j=0; j<lfsu_s.size(); j++)
959 for (size_type i=0; i<lfsu_s.size(); i++)
960 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * normalflux * factor * phi_s[i]);
967 std::vector<JacobianType> gradphi_s(lfsu_s.size());
968 lfsu_s.finiteElement().localBasis().evaluateJacobian(iplocal_s,gradphi_s);
970 const std::vector<JacobianType>& gradphi_s = cache[order_s].evaluateJacobian(iplocal_s,lfsu_s.finiteElement().localBasis());
974 jac = ig.inside()->geometry().jacobianInverseTransposed(iplocal_s);
975 std::vector<Dune::FieldVector<RF,dim> > tgradphi_s(lfsu_s.size());
976 for (size_type i=0; i<lfsu_s.size(); i++) jac.mv(gradphi_s[i][0],tgradphi_s[i]);
979 RF omegaup_s, omegaup_n;
992 for (size_type j=0; j<lfsu_s.size(); j++)
993 for (size_type i=0; i<lfsu_s.size(); i++)
994 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,omegaup_s * phi_s[j] * normalflux * factor * phi_s[i]);
997 for (size_type j=0; j<lfsu_s.size(); j++)
998 for (size_type i=0; i<lfsu_s.size(); i++)
999 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,-(An_F_s*tgradphi_s[j]) * factor * phi_s[i]);
1002 for (size_type j=0; j<lfsu_s.size(); j++)
1003 for (size_type i=0; i<lfsu_s.size(); i++)
1004 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,phi_s[j] * factor * theta * (An_F_s*tgradphi_s[i]));
1007 for (size_type j=0; j<lfsu_s.size(); j++)
1008 for (size_type i=0; i<lfsu_s.size(); i++)
1009 mat_ss.accumulate(lfsu_s,i,lfsu_s,j,penalty_factor * phi_s[j] * phi_s[i] * factor);
1014 template<
typename EG,
typename LFSV,
typename R>
1018 typedef typename LFSV::Traits::FiniteElementType::
1019 Traits::LocalBasisType::Traits::DomainFieldType DF;
1020 typedef typename LFSV::Traits::FiniteElementType::
1021 Traits::LocalBasisType::Traits::RangeFieldType RF;
1022 typedef typename LFSV::Traits::FiniteElementType::
1023 Traits::LocalBasisType::Traits::RangeType RangeType;
1024 typedef typename LFSV::Traits::SizeType size_type;
1027 const int dim = EG::Geometry::dimension;
1028 const int order = lfsv.finiteElement().localBasis().order();
1029 const int intorder = intorderadd + 2 * order;
1032 Dune::GeometryType gt = eg.geometry().type();
1033 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
1036 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
1040 std::vector<RangeType> phi(lfsv.size());
1041 lfsv.finiteElement().localBasis().evaluateFunction(it->position(),phi);
1043 const std::vector<RangeType>& phi = cache[order].evaluateFunction(it->position(),lfsv.finiteElement().localBasis());
1048 f = param.f(eg.entity(),it->position());
1051 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
1052 for (size_type i=0; i<lfsv.size(); i++)
1053 r.accumulate(lfsv,i,-f*phi[i]*factor);
1069 int quadrature_factor;
1071 typedef typename FiniteElementMap::Traits::FiniteElementType::Traits::LocalBasisType LocalBasisType;
1084 std::vector<Cache> cache;
1087 void element_size (
const GEO& geo,
typename GEO::ctype& hmin,
typename GEO::ctype hmax)
const
1089 typedef typename GEO::ctype DF;
1092 const int dim = GEO::coorddimension;
1095 Dune::FieldVector<DF,dim> x = geo.corner(0);
1097 hmin = hmax = x.two_norm();
1102 Dune::GeometryType gt = geo.type();
1103 for (
int i=0; i<Dune::ReferenceElements<DF,dim>::general(gt).size(dim-1); i++)
1105 Dune::FieldVector<DF,dim> x = geo.corner(Dune::ReferenceElements<DF,dim>::general(gt).subEntity(i,dim-1,0,dim));
1106 x -= geo.corner(Dune::ReferenceElements<DF,dim>::general(gt).subEntity(i,dim-1,1,dim));
1107 hmin = std::min(hmin,x.two_norm());
1108 hmax = std::max(hmax,x.two_norm());
Definition: convectiondiffusionparameter.hh:68
Definition: convectiondiffusiondg.hh:35
sparsity pattern generator
Definition: pattern.hh:30
Implement jacobian_apply_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:445
Definition: convectiondiffusiondg.hh:74
Default flags for all local operators.
Definition: flags.hh:18
Definition: convectiondiffusiondg.hh:80
Type
Definition: convectiondiffusiondg.hh:40
Implement jacobian_apply_boundary() based on alpha_boundary()
Definition: defaultimp.hh:538
void setTime(double t)
set time in parameter class
Definition: convectiondiffusiondg.hh:1058
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:276
Definition: convectiondiffusiondg.hh:40
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: convectiondiffusiondg.hh:1015
Definition: convectiondiffusiondg.hh:78
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: convectiondiffusiondg.hh:103
Definition: convectiondiffusiondg.hh:58
const IG & ig
Definition: common/constraints.hh:146
Definition: convectiondiffusiondg.hh:40
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: convectiondiffusiondg.hh:198
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:859
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:667
Definition: convectiondiffusiondg.hh:38
store values of basis functions and gradients in a cache
Definition: localbasiscache.hh:16
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
const F & f
Definition: common/constraints.hh:145
Definition: convectiondiffusiondg.hh:75
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:481
sparsity pattern generator
Definition: pattern.hh:14
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
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:84
const EG & eg
Definition: common/constraints.hh:277
Definition: convectiondiffusiondg.hh:33
Definition: convectiondiffusiondg.hh:79
Type
Definition: convectiondiffusiondg.hh:35
Definition: convectiondiffusiondg.hh:81
Definition: convectiondiffusionparameter.hh:68
const E & e
Definition: interpolate.hh:172
Type
Definition: convectiondiffusionparameter.hh:68
Definition: convectiondiffusionparameter.hh:68
Definition: convectiondiffusiondg.hh:35