2 #ifndef DUNE_PDELAB_DIFFUSIONDG_HH
3 #define DUNE_PDELAB_DIFFUSIONDG_HH
5 #include <dune/common/exceptions.hh>
6 #include <dune/common/fvector.hh>
7 #include <dune/geometry/quadraturerules.hh>
8 #include <dune/geometry/referenceelements.hh>
31 template<
typename K,
typename F,
typename B,
typename G,
typename J>
37 #ifdef JacobianBasedAlphaX
42 #ifdef NumericalJacobianX
43 #ifdef JacobianBasedAlphaX
44 #error You have provide either the alpha_* or the jacobian_* methods...
64 DiffusionDG (
const K& k_,
const F& f_,
const B& bctype_,
const G& g_,
const J& j_,
int dg_method,
int _superintegration_order = 0) :
65 k(k_), f(f_), bctype(bctype_), g(g_), j(j_), superintegration_order(_superintegration_order)
76 else if (dg_method == 1)
80 beta = 2.0 - 0.5*F::Traits::dimDomain;
83 else if (dg_method == 2)
87 beta = 2.0 - 0.5*F::Traits::dimDomain;
91 #ifndef JacobianBasedAlphaX
93 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
94 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
97 typedef typename LFSU::Traits::FiniteElementType::
98 Traits::LocalBasisType::Traits::DomainFieldType DF;
99 typedef typename LFSU::Traits::FiniteElementType::
100 Traits::LocalBasisType::Traits::RangeFieldType RF;
101 typedef typename LFSU::Traits::FiniteElementType::
102 Traits::LocalBasisType::Traits::JacobianType JacobianType;
105 const int dim = EG::Geometry::dimension;
108 Dune::GeometryType gt = eg.geometry().type();
109 const int qorder = std::max ( 2 * ( (
int)lfsu.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
110 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
113 typename K::Traits::RangeType tensor(0.0);
114 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
115 k.evaluate(eg.entity(),localcenter,tensor);
118 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
121 std::vector<JacobianType> js(lfsu.size());
122 lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
125 const typename EG::Geometry::JacobianInverseTransposed jac =
126 eg.geometry().jacobianInverseTransposed(it->position());
127 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
128 for (
size_t i=0; i<lfsu.size(); i++)
131 jac.umv(js[i][0],gradphi[i]);
135 Dune::FieldVector<RF,dim> gradu(0.0);
136 for (
size_t i=0; i<lfsu.size(); i++)
137 gradu.axpy(x(lfsu,i),gradphi[i]);
140 Dune::FieldVector<RF,dim> Kgradu(0.0);
141 tensor.umv(gradu,Kgradu);
144 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
145 for (
size_t i=0; i<lfsu.size(); i++)
147 r.accumulate( lfsv, i, Kgradu*gradphi[i] * factor ) ;
154 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
156 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
157 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
158 R& r_s, R& r_n)
const
161 typedef typename LFSU::Traits::FiniteElementType::
162 Traits::LocalBasisType::Traits::DomainFieldType DF;
163 typedef typename LFSU::Traits::FiniteElementType::
164 Traits::LocalBasisType::Traits::RangeFieldType RF;
165 typedef typename LFSV::Traits::FiniteElementType::
166 Traits::LocalBasisType::Traits::RangeType RangeType;
167 typedef typename LFSV::Traits::FiniteElementType::
168 Traits::LocalBasisType::Traits::JacobianType JacobianType;
171 const int dim = IG::dimension;
172 const int dimw = IG::dimensionworld;
175 Dune::GeometryType gtface = ig.geometryInInside().type();
176 const int qorder = std::max( 0, std::max(
177 2 * ( (
int)lfsu_s.finiteElement().localBasis().order() - 1 ),
178 2 * ( (
int)lfsu_n.finiteElement().localBasis().order() - 1 ))) + superintegration_order;
179 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
182 const Dune::FieldVector<DF,IG::dimension-1>& face_center =
183 Dune::ReferenceElements<DF,IG::dimension-1>::
184 general(ig.geometry().type()).position(0,0);
185 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
188 const Dune::FieldVector<DF,IG::dimension>&
189 inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
190 const Dune::FieldVector<DF,IG::dimension>&
191 outside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.outside()->type()).position(0,0);
192 typename K::Traits::RangeType permeability_s(0.0);
193 typename K::Traits::RangeType permeability_n(0.0);
194 k.evaluate(*(ig.inside()),inside_local,permeability_s);
195 k.evaluate(*(ig.outside()),outside_local,permeability_n);
198 RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
201 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
204 Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(it->position());
205 Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(it->position());
208 std::vector<JacobianType> js_s(lfsv_s.size());
209 lfsv_s.finiteElement().localBasis().evaluateJacobian(local_s,js_s);
210 std::vector<JacobianType> js_n(lfsv_n.size());
211 lfsv_n.finiteElement().localBasis().evaluateJacobian(local_n,js_n);
214 typename IG::Entity::Geometry::JacobianInverseTransposed jac_s;
215 jac_s = ig.inside()->geometry().jacobianInverseTransposed(local_s);
216 std::vector<Dune::FieldVector<RF,dim> > gradphi_s(lfsv_s.size());
217 for (
size_t i=0; i<lfsv_s.size(); i++)
220 jac_s.umv(js_s[i][0],gradphi_s[i]);
222 typename IG::Entity::Geometry::JacobianInverseTransposed jac_n;
223 jac_n = ig.outside()->geometry().jacobianInverseTransposed(local_n);
224 std::vector<Dune::FieldVector<RF,dim> > gradphi_n(lfsv_n.size());
225 for (
size_t i=0; i<lfsv_n.size(); i++)
228 jac_n.umv(js_n[i][0],gradphi_n[i]);
232 std::vector<RangeType> phi_s(lfsv_s.size());
233 lfsv_s.finiteElement().localBasis().evaluateFunction(local_s,phi_s);
234 std::vector<RangeType> phi_n(lfsv_n.size());
235 lfsv_n.finiteElement().localBasis().evaluateFunction(local_n,phi_n);
238 Dune::FieldVector<RF,dim> gradu_s(0.0);
239 for (
size_t i=0; i<lfsu_s.size(); i++)
241 gradu_s.axpy(x_s(lfsu_s,i),gradphi_s[i]);
243 Dune::FieldVector<RF,dim> gradu_n(0.0);
244 for (
size_t i=0; i<lfsu_n.size(); i++)
246 gradu_n.axpy(x_n(lfsu_n,i),gradphi_n[i]);
250 Dune::FieldVector<RF,dim> kgradu_s(0.0);
251 permeability_s.umv(gradu_s,kgradu_s);
252 Dune::FieldVector<RF,dim> kgradu_n(0.0);
253 permeability_n.umv(gradu_n,kgradu_n);
257 for (
size_t i=0; i<lfsu_s.size(); i++)
259 u_s += x_s(lfsu_s,i)*phi_s[i];
262 for (
size_t i=0; i<lfsu_n.size(); i++)
264 u_n += x_n(lfsu_n,i)*phi_n[i];
268 RF u_jump = u_s - u_n;
271 RF kgradunormal_average = (kgradu_s + kgradu_n)*normal * 0.5;
274 std::vector<Dune::FieldVector<RF,dim> > kgradphi_s(lfsu_s.size());
275 std::vector<Dune::FieldVector<RF,dim> > kgradphi_n(lfsu_n.size());
276 for (
size_t i=0; i<lfsu_s.size(); i++)
278 permeability_s.mv(gradphi_s[i],kgradphi_s[i]);
280 for (
size_t i=0; i<lfsu_n.size(); i++)
282 permeability_n.mv(gradphi_n[i],kgradphi_n[i]);
286 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
287 for (
unsigned int i=0; i<lfsv_s.size(); i++)
290 r_s.accumulate( lfsv_s, i, penalty_weight * u_jump*phi_s[i]*factor );
292 r_s.accumulate( lfsv_s, i, epsilon*(kgradphi_s[i]*normal)*0.5*u_jump*factor );
293 r_s.accumulate( lfsv_s, i, - phi_s[i]*kgradunormal_average*factor );
295 for (
unsigned int i=0; i<lfsv_n.size(); i++)
298 r_n.accumulate( lfsv_s, i, penalty_weight * u_jump*(-phi_n[i])*factor );
300 r_n.accumulate( lfsv_s, i, epsilon*(kgradphi_n[i]*normal)*0.5*u_jump*factor );
301 r_n.accumulate( lfsv_s, i, phi_n[i] * kgradunormal_average * factor );
307 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
308 void alpha_boundary (
const IG&
ig,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
311 typedef typename LFSU::Traits::FiniteElementType::
312 Traits::LocalBasisType::Traits::DomainFieldType DF;
313 typedef typename LFSU::Traits::FiniteElementType::
314 Traits::LocalBasisType::Traits::RangeFieldType RF;
315 typedef typename LFSV::Traits::FiniteElementType::
316 Traits::LocalBasisType::Traits::RangeType RangeType;
317 typedef typename LFSV::Traits::FiniteElementType::
318 Traits::LocalBasisType::Traits::JacobianType JacobianType;
321 const int dim = IG::dimension;
322 const int dimw = IG::dimensionworld;
325 Dune::GeometryType gtface = ig.geometryInInside().type();
326 const int qorder = std::max ( 2 * ( (
int)lfsu.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
327 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
331 if( bctype.isDirichlet( ig, rule.begin()->position() ) )
334 const Dune::FieldVector<DF,IG::dimension-1>& face_center =
335 Dune::ReferenceElements<DF,IG::dimension-1>::
336 general(ig.geometry().type()).position(0,0);
338 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
341 const Dune::FieldVector<DF,IG::dimension>
342 localcenter = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
343 typename K::Traits::RangeType tensor(0.0);
344 k.evaluate(*ig.inside(),localcenter,tensor);
347 RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
350 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
353 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
356 std::vector<JacobianType> js(lfsv.size());
357 lfsv.finiteElement().localBasis().evaluateJacobian(local,js);
360 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
361 jac = ig.inside()->geometry().jacobianInverseTransposed(local);
362 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsv.size());
363 for (
size_t i=0; i<lfsv.size(); i++)
366 jac.umv(js[i][0],gradphi[i]);
370 std::vector<RangeType> phi(lfsv.size());
371 lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
374 Dune::FieldVector<RF,dim> gradu(0.0);
375 for (
size_t i=0; i<lfsu.size(); i++)
377 gradu.axpy(x(lfsu,i),gradphi[i]);
381 Dune::FieldVector<RF,dim> Kgradu(0.0);
382 tensor.umv(gradu,Kgradu);
386 for (
size_t i=0; i<lfsu.size(); i++)
388 u += x(lfsu,i)*phi[i];
392 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
393 for (
size_t i=0; i<lfsv.size(); i++)
396 r.accumulate( lfsv, i, penalty_weight*u*phi[i]*factor );
398 Dune::FieldVector<RF,dim> Kgradv(0.0);
399 tensor.umv(gradphi[i],Kgradv);
400 r.accumulate( lfsv, i, epsilon * (Kgradv*normal)*u*factor );
401 r.accumulate( lfsv, i, - phi[i]*(Kgradu*normal)*factor );
410 template<
typename EG,
typename LFSV,
typename R>
414 typedef typename LFSV::Traits::FiniteElementType::
415 Traits::LocalBasisType::Traits::DomainFieldType DF;
416 typedef typename LFSV::Traits::FiniteElementType::
417 Traits::LocalBasisType::Traits::RangeFieldType RF;
418 typedef typename LFSV::Traits::FiniteElementType::
419 Traits::LocalBasisType::Traits::RangeType RangeType;
422 const int dim = EG::Geometry::dimension;
425 Dune::GeometryType gt = eg.geometry().type();
426 const int qorder = std::max ( 2 * ( (
int)lfsv.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
427 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
430 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
433 std::vector<RangeType> phi(lfsv.size());
434 lfsv.finiteElement().localBasis().evaluateFunction(it->position(),phi);
437 typename F::Traits::RangeType y;
438 f.evaluate(eg.entity(),it->position(),y);
441 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
442 for (
size_t i=0; i<lfsv.size(); i++)
445 r.accumulate( lfsv, i, - y*phi[i]*factor );
452 template<
typename IG,
typename LFSV,
typename R>
456 typedef typename LFSV::Traits::FiniteElementType::
457 Traits::LocalBasisType::Traits::DomainFieldType DF;
458 typedef typename LFSV::Traits::FiniteElementType::
459 Traits::LocalBasisType::Traits::RangeFieldType RF;
460 typedef typename LFSV::Traits::FiniteElementType::
461 Traits::LocalBasisType::Traits::RangeType RangeType;
462 typedef typename LFSV::Traits::FiniteElementType::
463 Traits::LocalBasisType::Traits::JacobianType JacobianType;
466 const int dim = IG::dimension;
467 const int dimw = IG::dimensionworld;
470 Dune::GeometryType gtface = ig.geometryInInside().type();
471 const int qorder = std::max ( 2 * ( (
int)lfsv.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
472 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
476 if( bctype.isNeumann( ig, rule.begin()->position() ) )
479 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
482 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
485 std::vector<RangeType> phi(lfsv.size());
486 lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
489 typename J::Traits::RangeType y(0.0);
490 j.evaluate(*(ig.inside()),local,y);
493 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
494 for (
size_t i=0; i<lfsv.size(); i++)
497 r.accumulate( lfsv, i, y*phi[i]*factor );
502 else if( bctype.isDirichlet( ig, rule.begin()->position() ) )
508 const Dune::FieldVector<DF,IG::dimension-1>& face_center =
509 Dune::ReferenceElements<DF,IG::dimension-1>::
510 general(ig.geometry().type()).position(0,0);
512 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
514 const Dune::FieldVector<DF,IG::dimension>
515 localcenter = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
516 typename K::Traits::RangeType tensor(0.0);
517 k.evaluate(*ig.inside(),localcenter,tensor);
519 RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
522 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
525 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
528 std::vector<JacobianType> js(lfsv.size());
529 lfsv.finiteElement().localBasis().evaluateJacobian(local,js);
532 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
533 jac = ig.inside()->geometry().jacobianInverseTransposed(local);
534 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsv.size());
535 for (
size_t i=0; i<lfsv.size(); i++)
538 jac.umv(js[i][0],gradphi[i]);
542 std::vector<RangeType> phi(lfsv.size());
543 lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
546 typename G::Traits::RangeType y = 0;
547 g.evaluate(*(ig.inside()),local,y);
550 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
551 for (
size_t i=0; i<lfsv.size(); i++)
554 r.accumulate( lfsv, i, -penalty_weight * y * phi[i] * factor );
556 Dune::FieldVector<RF,dim> Kgradv(0.0);
557 tensor.umv(gradphi[i],Kgradv);
558 r.accumulate( lfsv, i, -epsilon * (Kgradv*normal)*y*factor );
564 DUNE_THROW( Dune::Exception,
"Unrecognized or unsupported boundary condition type!" );
568 #ifndef NumericalJacobianX
570 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
575 typedef typename LFSU::Traits::FiniteElementType::
576 Traits::LocalBasisType::Traits::DomainFieldType DF;
577 typedef typename LFSU::Traits::FiniteElementType::
578 Traits::LocalBasisType::Traits::RangeFieldType RF;
579 typedef typename LFSU::Traits::FiniteElementType::
580 Traits::LocalBasisType::Traits::JacobianType JacobianType;
583 const int dim = EG::Geometry::dimension;
586 Dune::GeometryType gt = eg.geometry().type();
587 const int qorder = std::max ( 2 * ( (
int)lfsu.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
588 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
591 typename K::Traits::RangeType tensor;
592 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
593 k.evaluate(eg.entity(),localcenter,tensor);
596 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
599 std::vector<JacobianType> js(lfsu.size());
600 lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
603 typename EG::Geometry::JacobianInverseTransposed jac;
604 jac = eg.geometry().jacobianInverseTransposed(it->position());
605 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
606 for (
typename LFSU::Traits::SizeType i=0; i<lfsu.size(); i++)
609 jac.umv(js[i][0],gradphi[i]);
613 std::vector<Dune::FieldVector<RF,dim> > Kgradphi(lfsu.size());
614 for (
typename LFSU::Traits::SizeType i=0; i<lfsu.size(); i++)
616 tensor.mv(gradphi[i],Kgradphi[i]);
620 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
621 for (
typename LFSU::Traits::SizeType j=0; j<lfsu.size(); j++)
623 for (
typename LFSU::Traits::SizeType i=0; i<lfsu.size(); i++)
626 mat.accumulate( lfsu, i, lfsu, j, ( Kgradphi[j]*gradphi[i])*factor );
633 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
635 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
636 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
637 M& mat_ss, M& mat_sn,
638 M& mat_ns, M& mat_nn)
const
641 typedef typename LFSU::Traits::FiniteElementType::
642 Traits::LocalBasisType::Traits::DomainFieldType DF;
643 typedef typename LFSU::Traits::FiniteElementType::
644 Traits::LocalBasisType::Traits::RangeFieldType RF;
645 typedef typename LFSV::Traits::FiniteElementType::
646 Traits::LocalBasisType::Traits::RangeType RangeType;
647 typedef typename LFSV::Traits::FiniteElementType::
648 Traits::LocalBasisType::Traits::JacobianType JacobianType;
651 const int dim = IG::dimension;
652 const int dimw = IG::dimensionworld;
655 Dune::GeometryType gtface = ig.geometryInInside().type();
656 const int qorder = std::max( 0, std::max(
657 2 * ( (
int)lfsu_s.finiteElement().localBasis().order() - 1 ),
658 2 * ( (
int)lfsu_n.finiteElement().localBasis().order() - 1 ))) + superintegration_order;
659 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
662 const Dune::FieldVector<DF,IG::dimension-1>& face_center =
663 Dune::ReferenceElements<DF,IG::dimension-1>::
664 general(ig.geometry().type()).position(0,0);
665 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
668 const Dune::FieldVector<DF,IG::dimension>&
669 inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
670 const Dune::FieldVector<DF,IG::dimension>&
671 outside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.outside()->type()).position(0,0);
672 typename K::Traits::RangeType permeability_s(0.0);
673 typename K::Traits::RangeType permeability_n(0.0);
674 k.evaluate(*(ig.inside()),inside_local,permeability_s);
675 k.evaluate(*(ig.outside()),outside_local,permeability_n);
678 RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
681 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
684 Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(it->position());
685 Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(it->position());
688 std::vector<JacobianType> js_s(lfsv_s.size());
689 lfsv_s.finiteElement().localBasis().evaluateJacobian(local_s,js_s);
690 std::vector<JacobianType> js_n(lfsv_n.size());
691 lfsv_n.finiteElement().localBasis().evaluateJacobian(local_n,js_n);
694 typename IG::Entity::Geometry::JacobianInverseTransposed jac_s;
695 jac_s = ig.inside()->geometry().jacobianInverseTransposed(local_s);
696 std::vector<Dune::FieldVector<RF,dim> > gradphi_s(lfsv_s.size());
697 for (
size_t i=0; i<lfsv_s.size(); i++)
700 jac_s.umv(js_s[i][0],gradphi_s[i]);
702 typename IG::Entity::Geometry::JacobianInverseTransposed jac_n;
703 jac_n = ig.outside()->geometry().jacobianInverseTransposed(local_n);
704 std::vector<Dune::FieldVector<RF,dim> > gradphi_n(lfsv_n.size());
705 for (
size_t i=0; i<lfsv_n.size(); i++)
708 jac_n.umv(js_n[i][0],gradphi_n[i]);
712 std::vector<RangeType> phi_s(lfsv_s.size());
713 lfsv_s.finiteElement().localBasis().evaluateFunction(local_s,phi_s);
714 std::vector<RangeType> phi_n(lfsv_n.size());
715 lfsv_n.finiteElement().localBasis().evaluateFunction(local_n,phi_n);
718 Dune::FieldVector<RF,dim> gradu_s(0.0);
719 for (
size_t i=0; i<lfsu_s.size(); i++)
721 gradu_s.axpy(x_s(lfsu_s,i),gradphi_s[i]);
723 Dune::FieldVector<RF,dim> gradu_n(0.0);
724 for (
size_t i=0; i<lfsu_n.size(); i++)
726 gradu_n.axpy(x_n(lfsu_n,i),gradphi_n[i]);
730 Dune::FieldVector<RF,dim> kgradu_s(0.0);
731 permeability_s.umv(gradu_s,kgradu_s);
732 Dune::FieldVector<RF,dim> kgradu_n(0.0);
733 permeability_n.umv(gradu_n,kgradu_n);
737 for (
size_t i=0; i<lfsu_s.size(); i++)
739 u_s += x_s(lfsu_n,i)*phi_s[i];
742 for (
size_t i=0; i<lfsu_n.size(); i++)
744 u_n += x_n(lfsu_n,i)*phi_n[i];
748 std::vector<Dune::FieldVector<RF,dim> > kgradphi_s(lfsu_s.size());
749 std::vector<Dune::FieldVector<RF,dim> > kgradphi_n(lfsu_n.size());
750 for (
size_t i=0; i<lfsu_s.size(); i++)
752 permeability_s.mv(gradphi_s[i],kgradphi_s[i]);
754 for (
size_t i=0; i<lfsu_n.size(); i++)
756 permeability_n.mv(gradphi_n[i],kgradphi_n[i]);
760 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
761 for (
typename LFSU::Traits::SizeType j=0; j<lfsu_s.size(); j++)
763 for (
typename LFSU::Traits::SizeType i=0; i<lfsu_s.size(); i++)
766 mat_ss.accumulate( lfsu_s, i, lfsu_s, j, (epsilon*0.5*(kgradphi_s[i]*normal)*(phi_s[j]) - (phi_s[i])*0.5*(kgradphi_s[j]*normal) ) * factor );
768 mat_ss.accumulate( lfsu_s, i, lfsu_s, j, penalty_weight*(phi_s[j])*(phi_s[i]) * factor );
770 for (
typename LFSU::Traits::SizeType i=0; i<lfsu_n.size(); i++)
773 mat_ns.accumulate( lfsu_n, i, lfsu_s, j, (epsilon*0.5*(kgradphi_n[i]*normal)*(phi_s[j]) - (-phi_n[i])*0.5*(kgradphi_s[j]*normal) )*factor );
775 mat_ns.accumulate( lfsu_n, i, lfsu_s, j, penalty_weight*(phi_s[j])*(-phi_n[i]) *factor );
778 for (
typename LFSU::Traits::SizeType j=0; j<lfsu_n.size(); j++)
780 for (
typename LFSU::Traits::SizeType i=0; i<lfsu_s.size(); i++)
783 mat_sn.accumulate( lfsu_s, i, lfsu_n, j, (epsilon*0.5*(kgradphi_s[i]*normal)*(-phi_n[j]) - (phi_s[i])*0.5*(kgradphi_n[j]*normal) )*factor );
785 mat_sn.accumulate( lfsu_s, i, lfsu_n, j, penalty_weight*(-phi_n[j])*(phi_s[i]) *factor );
787 for (
typename LFSU::Traits::SizeType i=0; i<lfsu_n.size(); i++)
790 mat_nn.accumulate( lfsu_s, i, lfsu_n, j, (epsilon*0.5*(kgradphi_n[i]*normal)*(-phi_n[j]) - (-phi_n[i])*0.5*(kgradphi_n[j]*normal) )*factor );
792 mat_nn.accumulate( lfsu_s, i, lfsu_n, j, penalty_weight*(-phi_n[j])*(-phi_n[i]) *factor );
799 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
801 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
805 typedef typename LFSU::Traits::FiniteElementType::
806 Traits::LocalBasisType::Traits::DomainFieldType DF;
807 typedef typename LFSU::Traits::FiniteElementType::
808 Traits::LocalBasisType::Traits::RangeFieldType RF;
809 typedef typename LFSV::Traits::FiniteElementType::
810 Traits::LocalBasisType::Traits::RangeType RangeType;
811 typedef typename LFSV::Traits::FiniteElementType::
812 Traits::LocalBasisType::Traits::JacobianType JacobianType;
815 const int dim = IG::dimension;
816 const int dimw = IG::dimensionworld;
819 Dune::GeometryType gtface = ig.geometryInInside().type();
820 const int qorder = std::max ( 2 * ( (
int)lfsu.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
821 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
825 if( bctype.isDirichlet( ig, rule.begin()->position() ) )
828 const Dune::FieldVector<DF,IG::dimension-1>& face_center =
829 Dune::ReferenceElements<DF,IG::dimension-1>::
830 general(ig.geometry().type()).position(0,0);
832 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
835 const Dune::FieldVector<DF,IG::dimension>
836 localcenter = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
837 typename K::Traits::RangeType tensor(0.0);
838 k.evaluate(*ig.inside(),localcenter,tensor);
841 RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
844 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
847 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
850 std::vector<JacobianType> js(lfsv.size());
851 lfsv.finiteElement().localBasis().evaluateJacobian(local,js);
854 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
855 jac = ig.inside()->geometry().jacobianInverseTransposed(local);
856 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsv.size());
857 for (
size_t i=0; i<lfsv.size(); i++)
860 jac.umv(js[i][0],gradphi[i]);
864 std::vector<RangeType> phi(lfsv.size());
865 lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
868 Dune::FieldVector<RF,dim> gradu(0.0);
869 for (
size_t i=0; i<lfsu.size(); i++)
871 gradu.axpy(x(lfsu,i),gradphi[i]);
875 std::vector<Dune::FieldVector<RF,dim> > kgradphi(lfsu.size());
876 for (
size_t i=0; i<lfsu.size(); i++)
878 tensor.mv(gradphi[i],kgradphi[i]);
882 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
883 for (
typename LFSU::Traits::SizeType j=0; j<lfsu.size(); j++)
885 for (
typename LFSU::Traits::SizeType i=0; i<lfsu.size(); i++)
888 mat.accumulate( lfsu, i, lfsu, j, (epsilon*(kgradphi[i]*normal)*phi[j] - phi[i]*(kgradphi[j]*normal))*factor );
890 mat.accumulate( lfsu, i, lfsu, j, (penalty_weight*phi[j]*phi[i])*factor );
908 int superintegration_order;
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: diffusiondg.hh:634
Definition: diffusiondg.hh:32
void alpha_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusiondg.hh:308
Implement alpha_volume() based on jacobian_volume()
Definition: defaultimp.hh:611
Implement alpha_boundary() based on jacobian_boundary()
Definition: defaultimp.hh:704
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: diffusiondg.hh:155
void jacobian_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: diffusiondg.hh:800
DiffusionDG(const K &k_, const F &f_, const B &bctype_, const G &g_, const J &j_, int dg_method, int _superintegration_order=0)
Definition: diffusiondg.hh:64
Definition: diffusiondg.hh:62
Definition: diffusiondg.hh:54
const IG & ig
Definition: constraints.hh:147
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: diffusiondg.hh:453
Definition: diffusiondg.hh:59
Definition: diffusiondg.hh:57
Default flags for all local operators.
Definition: flags.hh:18
sparsity pattern generator
Definition: pattern.hh:13
Definition: diffusiondg.hh:60
Definition: diffusiondg.hh:61
Definition: diffusiondg.hh:53
Definition: adaptivity.hh:27
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: diffusiondg.hh:571
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
Implement jacobian_boundary() based on alpha_boundary()
Definition: defaultimp.hh:251
static const int dim
Definition: adaptivity.hh:83
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: diffusiondg.hh:411
sparsity pattern generator
Definition: pattern.hh:29
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusiondg.hh:94
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:157
Definition: diffusiondg.hh:58
Implement alpha_skeleton() based on jacobian_skeleton()
Definition: defaultimp.hh:650
const EG & eg
Definition: constraints.hh:280