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/common/static_assert.hh>
8 #include <dune/geometry/quadraturerules.hh>
9 #include <dune/geometry/referenceelements.hh>
32 template<
typename K,
typename F,
typename B,
typename G,
typename J>
38 #ifdef JacobianBasedAlphaX
43 #ifdef NumericalJacobianX
44 #ifdef JacobianBasedAlphaX
45 #error You have provide either the alpha_* or the jacobian_* methods...
65 DiffusionDG (
const K& k_,
const F& f_,
const B& bctype_,
const G& g_,
const J& j_,
int dg_method,
int _superintegration_order = 0) :
66 k(k_), f(f_), bctype(bctype_), g(g_), j(j_), superintegration_order(_superintegration_order)
77 else if (dg_method == 1)
81 beta = 2.0 - 0.5*F::Traits::dimDomain;
84 else if (dg_method == 2)
88 beta = 2.0 - 0.5*F::Traits::dimDomain;
92 #ifndef JacobianBasedAlphaX
94 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
95 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
98 typedef typename LFSU::Traits::FiniteElementType::
99 Traits::LocalBasisType::Traits::DomainFieldType DF;
100 typedef typename LFSU::Traits::FiniteElementType::
101 Traits::LocalBasisType::Traits::RangeFieldType RF;
102 typedef typename LFSU::Traits::FiniteElementType::
103 Traits::LocalBasisType::Traits::JacobianType JacobianType;
106 const int dim = EG::Geometry::dimension;
109 Dune::GeometryType gt = eg.geometry().type();
110 const int qorder = std::max ( 2 * ( (
int)lfsu.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
111 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
114 typename K::Traits::RangeType tensor(0.0);
115 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
116 k.evaluate(eg.entity(),localcenter,tensor);
119 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
122 std::vector<JacobianType> js(lfsu.size());
123 lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
126 const typename EG::Geometry::JacobianInverseTransposed jac =
127 eg.geometry().jacobianInverseTransposed(it->position());
128 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
129 for (
size_t i=0; i<lfsu.size(); i++)
132 jac.umv(js[i][0],gradphi[i]);
136 Dune::FieldVector<RF,dim> gradu(0.0);
137 for (
size_t i=0; i<lfsu.size(); i++)
138 gradu.axpy(x(lfsu,i),gradphi[i]);
141 Dune::FieldVector<RF,dim> Kgradu(0.0);
142 tensor.umv(gradu,Kgradu);
145 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
146 for (
size_t i=0; i<lfsu.size(); i++)
148 r.accumulate( lfsv, i, Kgradu*gradphi[i] * factor ) ;
155 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
157 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
158 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
159 R& r_s, R& r_n)
const
162 typedef typename LFSU::Traits::FiniteElementType::
163 Traits::LocalBasisType::Traits::DomainFieldType DF;
164 typedef typename LFSU::Traits::FiniteElementType::
165 Traits::LocalBasisType::Traits::RangeFieldType RF;
166 typedef typename LFSV::Traits::FiniteElementType::
167 Traits::LocalBasisType::Traits::RangeType RangeType;
168 typedef typename LFSV::Traits::FiniteElementType::
169 Traits::LocalBasisType::Traits::JacobianType JacobianType;
172 const int dim = IG::dimension;
173 const int dimw = IG::dimensionworld;
176 Dune::GeometryType gtface = ig.geometryInInside().type();
177 const int qorder = std::max( 0, std::max(
178 2 * ( (
int)lfsu_s.finiteElement().localBasis().order() - 1 ),
179 2 * ( (
int)lfsu_n.finiteElement().localBasis().order() - 1 ))) + superintegration_order;
180 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
183 const Dune::FieldVector<DF,IG::dimension-1>& face_center =
184 Dune::ReferenceElements<DF,IG::dimension-1>::
185 general(ig.geometry().type()).position(0,0);
186 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
189 const Dune::FieldVector<DF,IG::dimension>&
190 inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
191 const Dune::FieldVector<DF,IG::dimension>&
192 outside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.outside()->type()).position(0,0);
193 typename K::Traits::RangeType permeability_s(0.0);
194 typename K::Traits::RangeType permeability_n(0.0);
195 k.evaluate(*(ig.inside()),inside_local,permeability_s);
196 k.evaluate(*(ig.outside()),outside_local,permeability_n);
199 RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
202 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
205 Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(it->position());
206 Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(it->position());
209 std::vector<JacobianType> js_s(lfsv_s.size());
210 lfsv_s.finiteElement().localBasis().evaluateJacobian(local_s,js_s);
211 std::vector<JacobianType> js_n(lfsv_n.size());
212 lfsv_n.finiteElement().localBasis().evaluateJacobian(local_n,js_n);
215 typename IG::Entity::Geometry::JacobianInverseTransposed jac_s;
216 jac_s = ig.inside()->geometry().jacobianInverseTransposed(local_s);
217 std::vector<Dune::FieldVector<RF,dim> > gradphi_s(lfsv_s.size());
218 for (
size_t i=0; i<lfsv_s.size(); i++)
221 jac_s.umv(js_s[i][0],gradphi_s[i]);
223 typename IG::Entity::Geometry::JacobianInverseTransposed jac_n;
224 jac_n = ig.outside()->geometry().jacobianInverseTransposed(local_n);
225 std::vector<Dune::FieldVector<RF,dim> > gradphi_n(lfsv_n.size());
226 for (
size_t i=0; i<lfsv_n.size(); i++)
229 jac_n.umv(js_n[i][0],gradphi_n[i]);
233 std::vector<RangeType> phi_s(lfsv_s.size());
234 lfsv_s.finiteElement().localBasis().evaluateFunction(local_s,phi_s);
235 std::vector<RangeType> phi_n(lfsv_n.size());
236 lfsv_n.finiteElement().localBasis().evaluateFunction(local_n,phi_n);
239 Dune::FieldVector<RF,dim> gradu_s(0.0);
240 for (
size_t i=0; i<lfsu_s.size(); i++)
242 gradu_s.axpy(x_s(lfsu_s,i),gradphi_s[i]);
244 Dune::FieldVector<RF,dim> gradu_n(0.0);
245 for (
size_t i=0; i<lfsu_n.size(); i++)
247 gradu_n.axpy(x_n(lfsu_n,i),gradphi_n[i]);
251 Dune::FieldVector<RF,dim> kgradu_s(0.0);
252 permeability_s.umv(gradu_s,kgradu_s);
253 Dune::FieldVector<RF,dim> kgradu_n(0.0);
254 permeability_n.umv(gradu_n,kgradu_n);
258 for (
size_t i=0; i<lfsu_s.size(); i++)
260 u_s += x_s(lfsu_s,i)*phi_s[i];
263 for (
size_t i=0; i<lfsu_n.size(); i++)
265 u_n += x_n(lfsu_n,i)*phi_n[i];
269 RF u_jump = u_s - u_n;
272 RF kgradunormal_average = (kgradu_s + kgradu_n)*normal * 0.5;
275 std::vector<Dune::FieldVector<RF,dim> > kgradphi_s(lfsu_s.size());
276 std::vector<Dune::FieldVector<RF,dim> > kgradphi_n(lfsu_n.size());
277 for (
size_t i=0; i<lfsu_s.size(); i++)
279 permeability_s.mv(gradphi_s[i],kgradphi_s[i]);
281 for (
size_t i=0; i<lfsu_n.size(); i++)
283 permeability_n.mv(gradphi_n[i],kgradphi_n[i]);
287 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
288 for (
unsigned int i=0; i<lfsv_s.size(); i++)
291 r_s.accumulate( lfsv_s, i, penalty_weight * u_jump*phi_s[i]*factor );
293 r_s.accumulate( lfsv_s, i, epsilon*(kgradphi_s[i]*normal)*0.5*u_jump*factor );
294 r_s.accumulate( lfsv_s, i, - phi_s[i]*kgradunormal_average*factor );
296 for (
unsigned int i=0; i<lfsv_n.size(); i++)
299 r_n.accumulate( lfsv_s, i, penalty_weight * u_jump*(-phi_n[i])*factor );
301 r_n.accumulate( lfsv_s, i, epsilon*(kgradphi_n[i]*normal)*0.5*u_jump*factor );
302 r_n.accumulate( lfsv_s, i, phi_n[i] * kgradunormal_average * factor );
308 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
309 void alpha_boundary (
const IG&
ig,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
312 typedef typename LFSU::Traits::FiniteElementType::
313 Traits::LocalBasisType::Traits::DomainFieldType DF;
314 typedef typename LFSU::Traits::FiniteElementType::
315 Traits::LocalBasisType::Traits::RangeFieldType RF;
316 typedef typename LFSV::Traits::FiniteElementType::
317 Traits::LocalBasisType::Traits::RangeType RangeType;
318 typedef typename LFSV::Traits::FiniteElementType::
319 Traits::LocalBasisType::Traits::JacobianType JacobianType;
322 const int dim = IG::dimension;
323 const int dimw = IG::dimensionworld;
326 Dune::GeometryType gtface = ig.geometryInInside().type();
327 const int qorder = std::max ( 2 * ( (
int)lfsu.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
328 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
332 if( bctype.isDirichlet( ig, rule.begin()->position() ) )
335 const Dune::FieldVector<DF,IG::dimension-1>& face_center =
336 Dune::ReferenceElements<DF,IG::dimension-1>::
337 general(ig.geometry().type()).position(0,0);
339 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
342 const Dune::FieldVector<DF,IG::dimension>
343 localcenter = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
344 typename K::Traits::RangeType tensor(0.0);
345 k.evaluate(*ig.inside(),localcenter,tensor);
348 RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
351 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
354 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
357 std::vector<JacobianType> js(lfsv.size());
358 lfsv.finiteElement().localBasis().evaluateJacobian(local,js);
361 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
362 jac = ig.inside()->geometry().jacobianInverseTransposed(local);
363 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsv.size());
364 for (
size_t i=0; i<lfsv.size(); i++)
367 jac.umv(js[i][0],gradphi[i]);
371 std::vector<RangeType> phi(lfsv.size());
372 lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
375 Dune::FieldVector<RF,dim> gradu(0.0);
376 for (
size_t i=0; i<lfsu.size(); i++)
378 gradu.axpy(x(lfsu,i),gradphi[i]);
382 Dune::FieldVector<RF,dim> Kgradu(0.0);
383 tensor.umv(gradu,Kgradu);
387 for (
size_t i=0; i<lfsu.size(); i++)
389 u += x(lfsu,i)*phi[i];
393 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
394 for (
size_t i=0; i<lfsv.size(); i++)
397 r.accumulate( lfsv, i, penalty_weight*u*phi[i]*factor );
399 Dune::FieldVector<RF,dim> Kgradv(0.0);
400 tensor.umv(gradphi[i],Kgradv);
401 r.accumulate( lfsv, i, epsilon * (Kgradv*normal)*u*factor );
402 r.accumulate( lfsv, i, - phi[i]*(Kgradu*normal)*factor );
411 template<
typename EG,
typename LFSV,
typename R>
415 typedef typename LFSV::Traits::FiniteElementType::
416 Traits::LocalBasisType::Traits::DomainFieldType DF;
417 typedef typename LFSV::Traits::FiniteElementType::
418 Traits::LocalBasisType::Traits::RangeFieldType RF;
419 typedef typename LFSV::Traits::FiniteElementType::
420 Traits::LocalBasisType::Traits::RangeType RangeType;
423 const int dim = EG::Geometry::dimension;
426 Dune::GeometryType gt = eg.geometry().type();
427 const int qorder = std::max ( 2 * ( (
int)lfsv.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
428 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
431 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
434 std::vector<RangeType> phi(lfsv.size());
435 lfsv.finiteElement().localBasis().evaluateFunction(it->position(),phi);
438 typename F::Traits::RangeType y;
439 f.evaluate(eg.entity(),it->position(),y);
442 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
443 for (
size_t i=0; i<lfsv.size(); i++)
446 r.accumulate( lfsv, i, - y*phi[i]*factor );
453 template<
typename IG,
typename LFSV,
typename R>
457 typedef typename LFSV::Traits::FiniteElementType::
458 Traits::LocalBasisType::Traits::DomainFieldType DF;
459 typedef typename LFSV::Traits::FiniteElementType::
460 Traits::LocalBasisType::Traits::RangeFieldType RF;
461 typedef typename LFSV::Traits::FiniteElementType::
462 Traits::LocalBasisType::Traits::RangeType RangeType;
463 typedef typename LFSV::Traits::FiniteElementType::
464 Traits::LocalBasisType::Traits::JacobianType JacobianType;
467 const int dim = IG::dimension;
468 const int dimw = IG::dimensionworld;
471 Dune::GeometryType gtface = ig.geometryInInside().type();
472 const int qorder = std::max ( 2 * ( (
int)lfsv.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
473 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
477 if( bctype.isNeumann( ig, rule.begin()->position() ) )
480 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
483 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
486 std::vector<RangeType> phi(lfsv.size());
487 lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
490 typename J::Traits::RangeType y(0.0);
491 j.evaluate(*(ig.inside()),local,y);
494 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
495 for (
size_t i=0; i<lfsv.size(); i++)
498 r.accumulate( lfsv, i, y*phi[i]*factor );
503 else if( bctype.isDirichlet( ig, rule.begin()->position() ) )
509 const Dune::FieldVector<DF,IG::dimension-1>& face_center =
510 Dune::ReferenceElements<DF,IG::dimension-1>::
511 general(ig.geometry().type()).position(0,0);
513 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
515 const Dune::FieldVector<DF,IG::dimension>
516 localcenter = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
517 typename K::Traits::RangeType tensor(0.0);
518 k.evaluate(*ig.inside(),localcenter,tensor);
520 RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
523 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
526 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
529 std::vector<JacobianType> js(lfsv.size());
530 lfsv.finiteElement().localBasis().evaluateJacobian(local,js);
533 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
534 jac = ig.inside()->geometry().jacobianInverseTransposed(local);
535 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsv.size());
536 for (
size_t i=0; i<lfsv.size(); i++)
539 jac.umv(js[i][0],gradphi[i]);
543 std::vector<RangeType> phi(lfsv.size());
544 lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
547 typename G::Traits::RangeType y = 0;
548 g.evaluate(*(ig.inside()),local,y);
551 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
552 for (
size_t i=0; i<lfsv.size(); i++)
555 r.accumulate( lfsv, i, -penalty_weight * y * phi[i] * factor );
557 Dune::FieldVector<RF,dim> Kgradv(0.0);
558 tensor.umv(gradphi[i],Kgradv);
559 r.accumulate( lfsv, i, -epsilon * (Kgradv*normal)*y*factor );
565 DUNE_THROW( Dune::Exception,
"Unrecognized or unsupported boundary condition type!" );
569 #ifndef NumericalJacobianX
571 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
576 typedef typename LFSU::Traits::FiniteElementType::
577 Traits::LocalBasisType::Traits::DomainFieldType DF;
578 typedef typename LFSU::Traits::FiniteElementType::
579 Traits::LocalBasisType::Traits::RangeFieldType RF;
580 typedef typename LFSU::Traits::FiniteElementType::
581 Traits::LocalBasisType::Traits::JacobianType JacobianType;
584 const int dim = EG::Geometry::dimension;
587 Dune::GeometryType gt = eg.geometry().type();
588 const int qorder = std::max ( 2 * ( (
int)lfsu.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
589 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
592 typename K::Traits::RangeType tensor;
593 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
594 k.evaluate(eg.entity(),localcenter,tensor);
597 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
600 std::vector<JacobianType> js(lfsu.size());
601 lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
604 typename EG::Geometry::JacobianInverseTransposed jac;
605 jac = eg.geometry().jacobianInverseTransposed(it->position());
606 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
607 for (
typename LFSU::Traits::SizeType i=0; i<lfsu.size(); i++)
610 jac.umv(js[i][0],gradphi[i]);
614 std::vector<Dune::FieldVector<RF,dim> > Kgradphi(lfsu.size());
615 for (
typename LFSU::Traits::SizeType i=0; i<lfsu.size(); i++)
617 tensor.mv(gradphi[i],Kgradphi[i]);
621 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
622 for (
typename LFSU::Traits::SizeType j=0; j<lfsu.size(); j++)
624 for (
typename LFSU::Traits::SizeType i=0; i<lfsu.size(); i++)
627 mat.accumulate( lfsu, i, lfsu, j, ( Kgradphi[j]*gradphi[i])*factor );
634 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
636 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
637 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
638 M& mat_ss, M& mat_sn,
639 M& mat_ns, M& mat_nn)
const
642 typedef typename LFSU::Traits::FiniteElementType::
643 Traits::LocalBasisType::Traits::DomainFieldType DF;
644 typedef typename LFSU::Traits::FiniteElementType::
645 Traits::LocalBasisType::Traits::RangeFieldType RF;
646 typedef typename LFSV::Traits::FiniteElementType::
647 Traits::LocalBasisType::Traits::RangeType RangeType;
648 typedef typename LFSV::Traits::FiniteElementType::
649 Traits::LocalBasisType::Traits::JacobianType JacobianType;
652 const int dim = IG::dimension;
653 const int dimw = IG::dimensionworld;
656 Dune::GeometryType gtface = ig.geometryInInside().type();
657 const int qorder = std::max( 0, std::max(
658 2 * ( (
int)lfsu_s.finiteElement().localBasis().order() - 1 ),
659 2 * ( (
int)lfsu_n.finiteElement().localBasis().order() - 1 ))) + superintegration_order;
660 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
663 const Dune::FieldVector<DF,IG::dimension-1>& face_center =
664 Dune::ReferenceElements<DF,IG::dimension-1>::
665 general(ig.geometry().type()).position(0,0);
666 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
669 const Dune::FieldVector<DF,IG::dimension>&
670 inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
671 const Dune::FieldVector<DF,IG::dimension>&
672 outside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.outside()->type()).position(0,0);
673 typename K::Traits::RangeType permeability_s(0.0);
674 typename K::Traits::RangeType permeability_n(0.0);
675 k.evaluate(*(ig.inside()),inside_local,permeability_s);
676 k.evaluate(*(ig.outside()),outside_local,permeability_n);
679 RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
682 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
685 Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(it->position());
686 Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(it->position());
689 std::vector<JacobianType> js_s(lfsv_s.size());
690 lfsv_s.finiteElement().localBasis().evaluateJacobian(local_s,js_s);
691 std::vector<JacobianType> js_n(lfsv_n.size());
692 lfsv_n.finiteElement().localBasis().evaluateJacobian(local_n,js_n);
695 typename IG::Entity::Geometry::JacobianInverseTransposed jac_s;
696 jac_s = ig.inside()->geometry().jacobianInverseTransposed(local_s);
697 std::vector<Dune::FieldVector<RF,dim> > gradphi_s(lfsv_s.size());
698 for (
size_t i=0; i<lfsv_s.size(); i++)
701 jac_s.umv(js_s[i][0],gradphi_s[i]);
703 typename IG::Entity::Geometry::JacobianInverseTransposed jac_n;
704 jac_n = ig.outside()->geometry().jacobianInverseTransposed(local_n);
705 std::vector<Dune::FieldVector<RF,dim> > gradphi_n(lfsv_n.size());
706 for (
size_t i=0; i<lfsv_n.size(); i++)
709 jac_n.umv(js_n[i][0],gradphi_n[i]);
713 std::vector<RangeType> phi_s(lfsv_s.size());
714 lfsv_s.finiteElement().localBasis().evaluateFunction(local_s,phi_s);
715 std::vector<RangeType> phi_n(lfsv_n.size());
716 lfsv_n.finiteElement().localBasis().evaluateFunction(local_n,phi_n);
719 Dune::FieldVector<RF,dim> gradu_s(0.0);
720 for (
size_t i=0; i<lfsu_s.size(); i++)
722 gradu_s.axpy(x_s(lfsu_s,i),gradphi_s[i]);
724 Dune::FieldVector<RF,dim> gradu_n(0.0);
725 for (
size_t i=0; i<lfsu_n.size(); i++)
727 gradu_n.axpy(x_n(lfsu_n,i),gradphi_n[i]);
731 Dune::FieldVector<RF,dim> kgradu_s(0.0);
732 permeability_s.umv(gradu_s,kgradu_s);
733 Dune::FieldVector<RF,dim> kgradu_n(0.0);
734 permeability_n.umv(gradu_n,kgradu_n);
738 for (
size_t i=0; i<lfsu_s.size(); i++)
740 u_s += x_s(lfsu_n,i)*phi_s[i];
743 for (
size_t i=0; i<lfsu_n.size(); i++)
745 u_n += x_n(lfsu_n,i)*phi_n[i];
749 std::vector<Dune::FieldVector<RF,dim> > kgradphi_s(lfsu_s.size());
750 std::vector<Dune::FieldVector<RF,dim> > kgradphi_n(lfsu_n.size());
751 for (
size_t i=0; i<lfsu_s.size(); i++)
753 permeability_s.mv(gradphi_s[i],kgradphi_s[i]);
755 for (
size_t i=0; i<lfsu_n.size(); i++)
757 permeability_n.mv(gradphi_n[i],kgradphi_n[i]);
761 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
762 for (
typename LFSU::Traits::SizeType j=0; j<lfsu_s.size(); j++)
764 for (
typename LFSU::Traits::SizeType i=0; i<lfsu_s.size(); i++)
767 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 );
769 mat_ss.accumulate( lfsu_s, i, lfsu_s, j, penalty_weight*(phi_s[j])*(phi_s[i]) * factor );
771 for (
typename LFSU::Traits::SizeType i=0; i<lfsu_n.size(); i++)
774 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 );
776 mat_ns.accumulate( lfsu_n, i, lfsu_s, j, penalty_weight*(phi_s[j])*(-phi_n[i]) *factor );
779 for (
typename LFSU::Traits::SizeType j=0; j<lfsu_n.size(); j++)
781 for (
typename LFSU::Traits::SizeType i=0; i<lfsu_s.size(); i++)
784 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 );
786 mat_sn.accumulate( lfsu_s, i, lfsu_n, j, penalty_weight*(-phi_n[j])*(phi_s[i]) *factor );
788 for (
typename LFSU::Traits::SizeType i=0; i<lfsu_n.size(); i++)
791 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 );
793 mat_nn.accumulate( lfsu_s, i, lfsu_n, j, penalty_weight*(-phi_n[j])*(-phi_n[i]) *factor );
800 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
802 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
806 typedef typename LFSU::Traits::FiniteElementType::
807 Traits::LocalBasisType::Traits::DomainFieldType DF;
808 typedef typename LFSU::Traits::FiniteElementType::
809 Traits::LocalBasisType::Traits::RangeFieldType RF;
810 typedef typename LFSV::Traits::FiniteElementType::
811 Traits::LocalBasisType::Traits::RangeType RangeType;
812 typedef typename LFSV::Traits::FiniteElementType::
813 Traits::LocalBasisType::Traits::JacobianType JacobianType;
816 const int dim = IG::dimension;
817 const int dimw = IG::dimensionworld;
820 Dune::GeometryType gtface = ig.geometryInInside().type();
821 const int qorder = std::max ( 2 * ( (
int)lfsu.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
822 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
826 if( bctype.isDirichlet( ig, rule.begin()->position() ) )
829 const Dune::FieldVector<DF,IG::dimension-1>& face_center =
830 Dune::ReferenceElements<DF,IG::dimension-1>::
831 general(ig.geometry().type()).position(0,0);
833 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
836 const Dune::FieldVector<DF,IG::dimension>
837 localcenter = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
838 typename K::Traits::RangeType tensor(0.0);
839 k.evaluate(*ig.inside(),localcenter,tensor);
842 RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
845 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
848 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
851 std::vector<JacobianType> js(lfsv.size());
852 lfsv.finiteElement().localBasis().evaluateJacobian(local,js);
855 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
856 jac = ig.inside()->geometry().jacobianInverseTransposed(local);
857 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsv.size());
858 for (
size_t i=0; i<lfsv.size(); i++)
861 jac.umv(js[i][0],gradphi[i]);
865 std::vector<RangeType> phi(lfsv.size());
866 lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
869 Dune::FieldVector<RF,dim> gradu(0.0);
870 for (
size_t i=0; i<lfsu.size(); i++)
872 gradu.axpy(x(lfsu,i),gradphi[i]);
876 std::vector<Dune::FieldVector<RF,dim> > kgradphi(lfsu.size());
877 for (
size_t i=0; i<lfsu.size(); i++)
879 tensor.mv(gradphi[i],kgradphi[i]);
883 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
884 for (
typename LFSU::Traits::SizeType j=0; j<lfsu.size(); j++)
886 for (
typename LFSU::Traits::SizeType i=0; i<lfsu.size(); i++)
889 mat.accumulate( lfsu, i, lfsu, j, (epsilon*(kgradphi[i]*normal)*phi[j] - phi[i]*(kgradphi[j]*normal))*factor );
891 mat.accumulate( lfsu, i, lfsu, j, (penalty_weight*phi[j]*phi[i])*factor );
909 int superintegration_order;
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusiondg.hh:95
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: diffusiondg.hh:572
sparsity pattern generator
Definition: pattern.hh:30
Definition: diffusiondg.hh:54
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:156
Default flags for all local operators.
Definition: flags.hh:18
Implement jacobian_boundary() based on alpha_boundary()
Definition: defaultimp.hh:251
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:635
void alpha_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusiondg.hh:309
static const int dim
Definition: adaptivity.hh:82
Definition: diffusiondg.hh:61
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: diffusiondg.hh:412
Definition: diffusiondg.hh:55
Implement alpha_volume() based on jacobian_volume()
Definition: defaultimp.hh:611
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:65
Implement jacobian_volume() based on alpha_volume()
Definition: defaultimp.hh:32
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: diffusiondg.hh:454
Definition: diffusiondg.hh:59
Definition: diffusiondg.hh:60
const IG & ig
Definition: common/constraints.hh:146
Definition: diffusiondg.hh:58
Implement alpha_skeleton() based on jacobian_skeleton()
Definition: defaultimp.hh:650
Definition: diffusiondg.hh:62
sparsity pattern generator
Definition: pattern.hh:14
void jacobian_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: diffusiondg.hh:801
const EG & eg
Definition: common/constraints.hh:277
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: defaultimp.hh:157
Definition: diffusiondg.hh:33
Implement alpha_boundary() based on jacobian_boundary()
Definition: defaultimp.hh:704
Definition: diffusiondg.hh:63