2 #ifndef DUNE_PDELAB_DIFFUSIONDG_HH 3 #define DUNE_PDELAB_DIFFUSIONDG_HH 4 #warning This file is deprecated and will be removed after the Dune-PDELab 2.4 release! Use the ConvectionDiffusionDG local operator from dune/pdelab/localoperator/convectiondiffusiondg.hh instead! 6 #include <dune/common/exceptions.hh> 7 #include <dune/common/fvector.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 DUNE_DEPRECATED_MSG(
"Deprecated in Dune-PDELab 2.4, use the local operator ConvectionDiffusionDG instead!")
66 DiffusionDG (const K& k_, const F& f_, const B& bctype_, const G& g_, const J& j_,
int dg_method,
int _superintegration_order = 0) :
67 k(k_), f(f_), bctype(bctype_), g(g_), j(j_), superintegration_order(_superintegration_order)
78 else if (dg_method == 1)
82 beta = 2.0 - 0.5*F::Traits::dimDomain;
85 else if (dg_method == 2)
89 beta = 2.0 - 0.5*F::Traits::dimDomain;
93 #ifndef JacobianBasedAlphaX 95 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
96 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 99 typedef typename LFSU::Traits::FiniteElementType::
100 Traits::LocalBasisType::Traits::DomainFieldType DF;
101 typedef typename LFSU::Traits::FiniteElementType::
102 Traits::LocalBasisType::Traits::RangeFieldType RF;
103 typedef typename LFSU::Traits::FiniteElementType::
104 Traits::LocalBasisType::Traits::JacobianType JacobianType;
107 const int dim = EG::Geometry::dimension;
110 Dune::GeometryType gt = eg.geometry().type();
111 const int qorder = std::max ( 2 * ( (
int)lfsu.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
112 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
115 typename K::Traits::RangeType tensor(0.0);
116 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
117 k.evaluate(eg.entity(),localcenter,tensor);
120 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
123 std::vector<JacobianType> js(lfsu.size());
124 lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
127 const typename EG::Geometry::JacobianInverseTransposed jac =
128 eg.geometry().jacobianInverseTransposed(it->position());
129 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
130 for (
size_t i=0; i<lfsu.size(); i++)
133 jac.umv(js[i][0],gradphi[i]);
137 Dune::FieldVector<RF,dim> gradu(0.0);
138 for (
size_t i=0; i<lfsu.size(); i++)
139 gradu.axpy(x(lfsu,i),gradphi[i]);
142 Dune::FieldVector<RF,dim> Kgradu(0.0);
143 tensor.umv(gradu,Kgradu);
146 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
147 for (
size_t i=0; i<lfsu.size(); i++)
149 r.accumulate( lfsv, i, Kgradu*gradphi[i] * factor ) ;
156 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
158 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
159 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
160 R& r_s, R& r_n)
const 163 typedef typename LFSU::Traits::FiniteElementType::
164 Traits::LocalBasisType::Traits::DomainFieldType DF;
165 typedef typename LFSU::Traits::FiniteElementType::
166 Traits::LocalBasisType::Traits::RangeFieldType RF;
167 typedef typename LFSV::Traits::FiniteElementType::
168 Traits::LocalBasisType::Traits::RangeType RangeType;
169 typedef typename LFSV::Traits::FiniteElementType::
170 Traits::LocalBasisType::Traits::JacobianType JacobianType;
173 const int dim = IG::dimension;
174 const int dimw = IG::dimensionworld;
177 Dune::GeometryType gtface = ig.geometryInInside().type();
178 const int qorder = std::max( 0, std::max(
179 2 * ( (
int)lfsu_s.finiteElement().localBasis().order() - 1 ),
180 2 * ( (
int)lfsu_n.finiteElement().localBasis().order() - 1 ))) + superintegration_order;
181 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
184 const Dune::FieldVector<DF,IG::dimension-1>& face_center =
185 Dune::ReferenceElements<DF,IG::dimension-1>::
186 general(ig.geometry().type()).position(0,0);
187 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
190 const Dune::FieldVector<DF,IG::dimension>&
191 inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
192 const Dune::FieldVector<DF,IG::dimension>&
193 outside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.outside()->type()).position(0,0);
194 typename K::Traits::RangeType permeability_s(0.0);
195 typename K::Traits::RangeType permeability_n(0.0);
196 k.evaluate(*(ig.inside()),inside_local,permeability_s);
197 k.evaluate(*(ig.outside()),outside_local,permeability_n);
200 RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
203 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
206 Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(it->position());
207 Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(it->position());
210 std::vector<JacobianType> js_s(lfsv_s.size());
211 lfsv_s.finiteElement().localBasis().evaluateJacobian(local_s,js_s);
212 std::vector<JacobianType> js_n(lfsv_n.size());
213 lfsv_n.finiteElement().localBasis().evaluateJacobian(local_n,js_n);
216 typename IG::Entity::Geometry::JacobianInverseTransposed jac_s;
217 jac_s = ig.inside()->geometry().jacobianInverseTransposed(local_s);
218 std::vector<Dune::FieldVector<RF,dim> > gradphi_s(lfsv_s.size());
219 for (
size_t i=0; i<lfsv_s.size(); i++)
222 jac_s.umv(js_s[i][0],gradphi_s[i]);
224 typename IG::Entity::Geometry::JacobianInverseTransposed jac_n;
225 jac_n = ig.outside()->geometry().jacobianInverseTransposed(local_n);
226 std::vector<Dune::FieldVector<RF,dim> > gradphi_n(lfsv_n.size());
227 for (
size_t i=0; i<lfsv_n.size(); i++)
230 jac_n.umv(js_n[i][0],gradphi_n[i]);
234 std::vector<RangeType> phi_s(lfsv_s.size());
235 lfsv_s.finiteElement().localBasis().evaluateFunction(local_s,phi_s);
236 std::vector<RangeType> phi_n(lfsv_n.size());
237 lfsv_n.finiteElement().localBasis().evaluateFunction(local_n,phi_n);
240 Dune::FieldVector<RF,dim> gradu_s(0.0);
241 for (
size_t i=0; i<lfsu_s.size(); i++)
243 gradu_s.axpy(x_s(lfsu_s,i),gradphi_s[i]);
245 Dune::FieldVector<RF,dim> gradu_n(0.0);
246 for (
size_t i=0; i<lfsu_n.size(); i++)
248 gradu_n.axpy(x_n(lfsu_n,i),gradphi_n[i]);
252 Dune::FieldVector<RF,dim> kgradu_s(0.0);
253 permeability_s.umv(gradu_s,kgradu_s);
254 Dune::FieldVector<RF,dim> kgradu_n(0.0);
255 permeability_n.umv(gradu_n,kgradu_n);
259 for (
size_t i=0; i<lfsu_s.size(); i++)
261 u_s += x_s(lfsu_s,i)*phi_s[i];
264 for (
size_t i=0; i<lfsu_n.size(); i++)
266 u_n += x_n(lfsu_n,i)*phi_n[i];
270 RF u_jump = u_s - u_n;
273 RF kgradunormal_average = (kgradu_s + kgradu_n)*normal * 0.5;
276 std::vector<Dune::FieldVector<RF,dim> > kgradphi_s(lfsu_s.size());
277 std::vector<Dune::FieldVector<RF,dim> > kgradphi_n(lfsu_n.size());
278 for (
size_t i=0; i<lfsu_s.size(); i++)
280 permeability_s.mv(gradphi_s[i],kgradphi_s[i]);
282 for (
size_t i=0; i<lfsu_n.size(); i++)
284 permeability_n.mv(gradphi_n[i],kgradphi_n[i]);
288 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
289 for (
unsigned int i=0; i<lfsv_s.size(); i++)
292 r_s.accumulate( lfsv_s, i, penalty_weight * u_jump*phi_s[i]*factor );
294 r_s.accumulate( lfsv_s, i, epsilon*(kgradphi_s[i]*normal)*0.5*u_jump*factor );
295 r_s.accumulate( lfsv_s, i, - phi_s[i]*kgradunormal_average*factor );
297 for (
unsigned int i=0; i<lfsv_n.size(); i++)
300 r_n.accumulate( lfsv_s, i, penalty_weight * u_jump*(-phi_n[i])*factor );
302 r_n.accumulate( lfsv_s, i, epsilon*(kgradphi_n[i]*normal)*0.5*u_jump*factor );
303 r_n.accumulate( lfsv_s, i, phi_n[i] * kgradunormal_average * factor );
309 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename R>
310 void alpha_boundary (
const IG&
ig,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 313 typedef typename LFSU::Traits::FiniteElementType::
314 Traits::LocalBasisType::Traits::DomainFieldType DF;
315 typedef typename LFSU::Traits::FiniteElementType::
316 Traits::LocalBasisType::Traits::RangeFieldType RF;
317 typedef typename LFSV::Traits::FiniteElementType::
318 Traits::LocalBasisType::Traits::RangeType RangeType;
319 typedef typename LFSV::Traits::FiniteElementType::
320 Traits::LocalBasisType::Traits::JacobianType JacobianType;
323 const int dim = IG::dimension;
324 const int dimw = IG::dimensionworld;
327 Dune::GeometryType gtface = ig.geometryInInside().type();
328 const int qorder = std::max ( 2 * ( (
int)lfsu.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
329 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
333 if( bctype.isDirichlet( ig, rule.begin()->position() ) )
336 const Dune::FieldVector<DF,IG::dimension-1>& face_center =
337 Dune::ReferenceElements<DF,IG::dimension-1>::
338 general(ig.geometry().type()).position(0,0);
340 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
343 const Dune::FieldVector<DF,IG::dimension>
344 localcenter = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
345 typename K::Traits::RangeType tensor(0.0);
346 k.evaluate(*ig.inside(),localcenter,tensor);
349 RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
352 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
355 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
358 std::vector<JacobianType> js(lfsv.size());
359 lfsv.finiteElement().localBasis().evaluateJacobian(local,js);
362 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
363 jac = ig.inside()->geometry().jacobianInverseTransposed(local);
364 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsv.size());
365 for (
size_t i=0; i<lfsv.size(); i++)
368 jac.umv(js[i][0],gradphi[i]);
372 std::vector<RangeType> phi(lfsv.size());
373 lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
376 Dune::FieldVector<RF,dim> gradu(0.0);
377 for (
size_t i=0; i<lfsu.size(); i++)
379 gradu.axpy(x(lfsu,i),gradphi[i]);
383 Dune::FieldVector<RF,dim> Kgradu(0.0);
384 tensor.umv(gradu,Kgradu);
388 for (
size_t i=0; i<lfsu.size(); i++)
390 u += x(lfsu,i)*phi[i];
394 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
395 for (
size_t i=0; i<lfsv.size(); i++)
398 r.accumulate( lfsv, i, penalty_weight*u*phi[i]*factor );
400 Dune::FieldVector<RF,dim> Kgradv(0.0);
401 tensor.umv(gradphi[i],Kgradv);
402 r.accumulate( lfsv, i, epsilon * (Kgradv*normal)*u*factor );
403 r.accumulate( lfsv, i, - phi[i]*(Kgradu*normal)*factor );
412 template<
typename EG,
typename LFSV,
typename R>
416 typedef typename LFSV::Traits::FiniteElementType::
417 Traits::LocalBasisType::Traits::DomainFieldType DF;
418 typedef typename LFSV::Traits::FiniteElementType::
419 Traits::LocalBasisType::Traits::RangeFieldType RF;
420 typedef typename LFSV::Traits::FiniteElementType::
421 Traits::LocalBasisType::Traits::RangeType RangeType;
424 const int dim = EG::Geometry::dimension;
427 Dune::GeometryType gt = eg.geometry().type();
428 const int qorder = std::max ( 2 * ( (
int)lfsv.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
429 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
432 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
435 std::vector<RangeType> phi(lfsv.size());
436 lfsv.finiteElement().localBasis().evaluateFunction(it->position(),phi);
439 typename F::Traits::RangeType y;
440 f.evaluate(eg.entity(),it->position(),y);
443 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
444 for (
size_t i=0; i<lfsv.size(); i++)
447 r.accumulate( lfsv, i, - y*phi[i]*factor );
454 template<
typename IG,
typename LFSV,
typename R>
458 typedef typename LFSV::Traits::FiniteElementType::
459 Traits::LocalBasisType::Traits::DomainFieldType DF;
460 typedef typename LFSV::Traits::FiniteElementType::
461 Traits::LocalBasisType::Traits::RangeFieldType RF;
462 typedef typename LFSV::Traits::FiniteElementType::
463 Traits::LocalBasisType::Traits::RangeType RangeType;
464 typedef typename LFSV::Traits::FiniteElementType::
465 Traits::LocalBasisType::Traits::JacobianType JacobianType;
468 const int dim = IG::dimension;
469 const int dimw = IG::dimensionworld;
472 Dune::GeometryType gtface = ig.geometryInInside().type();
473 const int qorder = std::max ( 2 * ( (
int)lfsv.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
474 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
478 if( bctype.isNeumann( ig, rule.begin()->position() ) )
481 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
484 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
487 std::vector<RangeType> phi(lfsv.size());
488 lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
491 typename J::Traits::RangeType y(0.0);
492 j.evaluate(*(ig.inside()),local,y);
495 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
496 for (
size_t i=0; i<lfsv.size(); i++)
499 r.accumulate( lfsv, i, y*phi[i]*factor );
504 else if( bctype.isDirichlet( ig, rule.begin()->position() ) )
510 const Dune::FieldVector<DF,IG::dimension-1>& face_center =
511 Dune::ReferenceElements<DF,IG::dimension-1>::
512 general(ig.geometry().type()).position(0,0);
514 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
516 const Dune::FieldVector<DF,IG::dimension>
517 localcenter = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
518 typename K::Traits::RangeType tensor(0.0);
519 k.evaluate(*ig.inside(),localcenter,tensor);
521 RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
524 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
527 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
530 std::vector<JacobianType> js(lfsv.size());
531 lfsv.finiteElement().localBasis().evaluateJacobian(local,js);
534 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
535 jac = ig.inside()->geometry().jacobianInverseTransposed(local);
536 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsv.size());
537 for (
size_t i=0; i<lfsv.size(); i++)
540 jac.umv(js[i][0],gradphi[i]);
544 std::vector<RangeType> phi(lfsv.size());
545 lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
548 typename G::Traits::RangeType y = 0;
549 g.evaluate(*(ig.inside()),local,y);
552 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
553 for (
size_t i=0; i<lfsv.size(); i++)
556 r.accumulate( lfsv, i, -penalty_weight * y * phi[i] * factor );
558 Dune::FieldVector<RF,dim> Kgradv(0.0);
559 tensor.umv(gradphi[i],Kgradv);
560 r.accumulate( lfsv, i, -epsilon * (Kgradv*normal)*y*factor );
566 DUNE_THROW( Dune::Exception,
"Unrecognized or unsupported boundary condition type!" );
570 #ifndef NumericalJacobianX 572 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
577 typedef typename LFSU::Traits::FiniteElementType::
578 Traits::LocalBasisType::Traits::DomainFieldType DF;
579 typedef typename LFSU::Traits::FiniteElementType::
580 Traits::LocalBasisType::Traits::RangeFieldType RF;
581 typedef typename LFSU::Traits::FiniteElementType::
582 Traits::LocalBasisType::Traits::JacobianType JacobianType;
585 const int dim = EG::Geometry::dimension;
588 Dune::GeometryType gt = eg.geometry().type();
589 const int qorder = std::max ( 2 * ( (
int)lfsu.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
590 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
593 typename K::Traits::RangeType tensor;
594 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
595 k.evaluate(eg.entity(),localcenter,tensor);
598 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
601 std::vector<JacobianType> js(lfsu.size());
602 lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
605 typename EG::Geometry::JacobianInverseTransposed jac;
606 jac = eg.geometry().jacobianInverseTransposed(it->position());
607 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
608 for (
typename LFSU::Traits::SizeType i=0; i<lfsu.size(); i++)
611 jac.umv(js[i][0],gradphi[i]);
615 std::vector<Dune::FieldVector<RF,dim> > Kgradphi(lfsu.size());
616 for (
typename LFSU::Traits::SizeType i=0; i<lfsu.size(); i++)
618 tensor.mv(gradphi[i],Kgradphi[i]);
622 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
623 for (
typename LFSU::Traits::SizeType j=0; j<lfsu.size(); j++)
625 for (
typename LFSU::Traits::SizeType i=0; i<lfsu.size(); i++)
628 mat.accumulate( lfsu, i, lfsu, j, ( Kgradphi[j]*gradphi[i])*factor );
635 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
637 const LFSU& lfsu_s,
const X& x_s,
const LFSV& lfsv_s,
638 const LFSU& lfsu_n,
const X& x_n,
const LFSV& lfsv_n,
639 M& mat_ss, M& mat_sn,
640 M& mat_ns, M& mat_nn)
const 643 typedef typename LFSU::Traits::FiniteElementType::
644 Traits::LocalBasisType::Traits::DomainFieldType DF;
645 typedef typename LFSU::Traits::FiniteElementType::
646 Traits::LocalBasisType::Traits::RangeFieldType RF;
647 typedef typename LFSV::Traits::FiniteElementType::
648 Traits::LocalBasisType::Traits::RangeType RangeType;
649 typedef typename LFSV::Traits::FiniteElementType::
650 Traits::LocalBasisType::Traits::JacobianType JacobianType;
653 const int dim = IG::dimension;
654 const int dimw = IG::dimensionworld;
657 Dune::GeometryType gtface = ig.geometryInInside().type();
658 const int qorder = std::max( 0, std::max(
659 2 * ( (
int)lfsu_s.finiteElement().localBasis().order() - 1 ),
660 2 * ( (
int)lfsu_n.finiteElement().localBasis().order() - 1 ))) + superintegration_order;
661 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
664 const Dune::FieldVector<DF,IG::dimension-1>& face_center =
665 Dune::ReferenceElements<DF,IG::dimension-1>::
666 general(ig.geometry().type()).position(0,0);
667 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
670 const Dune::FieldVector<DF,IG::dimension>&
671 inside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
672 const Dune::FieldVector<DF,IG::dimension>&
673 outside_local = Dune::ReferenceElements<DF,IG::dimension>::general(ig.outside()->type()).position(0,0);
674 typename K::Traits::RangeType permeability_s(0.0);
675 typename K::Traits::RangeType permeability_n(0.0);
676 k.evaluate(*(ig.inside()),inside_local,permeability_s);
677 k.evaluate(*(ig.outside()),outside_local,permeability_n);
680 RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
683 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
686 Dune::FieldVector<DF,dim> local_s = ig.geometryInInside().global(it->position());
687 Dune::FieldVector<DF,dim> local_n = ig.geometryInOutside().global(it->position());
690 std::vector<JacobianType> js_s(lfsv_s.size());
691 lfsv_s.finiteElement().localBasis().evaluateJacobian(local_s,js_s);
692 std::vector<JacobianType> js_n(lfsv_n.size());
693 lfsv_n.finiteElement().localBasis().evaluateJacobian(local_n,js_n);
696 typename IG::Entity::Geometry::JacobianInverseTransposed jac_s;
697 jac_s = ig.inside()->geometry().jacobianInverseTransposed(local_s);
698 std::vector<Dune::FieldVector<RF,dim> > gradphi_s(lfsv_s.size());
699 for (
size_t i=0; i<lfsv_s.size(); i++)
702 jac_s.umv(js_s[i][0],gradphi_s[i]);
704 typename IG::Entity::Geometry::JacobianInverseTransposed jac_n;
705 jac_n = ig.outside()->geometry().jacobianInverseTransposed(local_n);
706 std::vector<Dune::FieldVector<RF,dim> > gradphi_n(lfsv_n.size());
707 for (
size_t i=0; i<lfsv_n.size(); i++)
710 jac_n.umv(js_n[i][0],gradphi_n[i]);
714 std::vector<RangeType> phi_s(lfsv_s.size());
715 lfsv_s.finiteElement().localBasis().evaluateFunction(local_s,phi_s);
716 std::vector<RangeType> phi_n(lfsv_n.size());
717 lfsv_n.finiteElement().localBasis().evaluateFunction(local_n,phi_n);
720 Dune::FieldVector<RF,dim> gradu_s(0.0);
721 for (
size_t i=0; i<lfsu_s.size(); i++)
723 gradu_s.axpy(x_s(lfsu_s,i),gradphi_s[i]);
725 Dune::FieldVector<RF,dim> gradu_n(0.0);
726 for (
size_t i=0; i<lfsu_n.size(); i++)
728 gradu_n.axpy(x_n(lfsu_n,i),gradphi_n[i]);
732 Dune::FieldVector<RF,dim> kgradu_s(0.0);
733 permeability_s.umv(gradu_s,kgradu_s);
734 Dune::FieldVector<RF,dim> kgradu_n(0.0);
735 permeability_n.umv(gradu_n,kgradu_n);
739 for (
size_t i=0; i<lfsu_s.size(); i++)
741 u_s += x_s(lfsu_n,i)*phi_s[i];
744 for (
size_t i=0; i<lfsu_n.size(); i++)
746 u_n += x_n(lfsu_n,i)*phi_n[i];
750 std::vector<Dune::FieldVector<RF,dim> > kgradphi_s(lfsu_s.size());
751 std::vector<Dune::FieldVector<RF,dim> > kgradphi_n(lfsu_n.size());
752 for (
size_t i=0; i<lfsu_s.size(); i++)
754 permeability_s.mv(gradphi_s[i],kgradphi_s[i]);
756 for (
size_t i=0; i<lfsu_n.size(); i++)
758 permeability_n.mv(gradphi_n[i],kgradphi_n[i]);
762 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
763 for (
typename LFSU::Traits::SizeType j=0; j<lfsu_s.size(); j++)
765 for (
typename LFSU::Traits::SizeType i=0; i<lfsu_s.size(); i++)
768 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 );
770 mat_ss.accumulate( lfsu_s, i, lfsu_s, j, penalty_weight*(phi_s[j])*(phi_s[i]) * factor );
772 for (
typename LFSU::Traits::SizeType i=0; i<lfsu_n.size(); i++)
775 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 );
777 mat_ns.accumulate( lfsu_n, i, lfsu_s, j, penalty_weight*(phi_s[j])*(-phi_n[i]) *factor );
780 for (
typename LFSU::Traits::SizeType j=0; j<lfsu_n.size(); j++)
782 for (
typename LFSU::Traits::SizeType i=0; i<lfsu_s.size(); i++)
785 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 );
787 mat_sn.accumulate( lfsu_s, i, lfsu_n, j, penalty_weight*(-phi_n[j])*(phi_s[i]) *factor );
789 for (
typename LFSU::Traits::SizeType i=0; i<lfsu_n.size(); i++)
792 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 );
794 mat_nn.accumulate( lfsu_s, i, lfsu_n, j, penalty_weight*(-phi_n[j])*(-phi_n[i]) *factor );
801 template<
typename IG,
typename LFSU,
typename X,
typename LFSV,
typename M>
803 const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
807 typedef typename LFSU::Traits::FiniteElementType::
808 Traits::LocalBasisType::Traits::DomainFieldType DF;
809 typedef typename LFSU::Traits::FiniteElementType::
810 Traits::LocalBasisType::Traits::RangeFieldType RF;
811 typedef typename LFSV::Traits::FiniteElementType::
812 Traits::LocalBasisType::Traits::RangeType RangeType;
813 typedef typename LFSV::Traits::FiniteElementType::
814 Traits::LocalBasisType::Traits::JacobianType JacobianType;
817 const int dim = IG::dimension;
818 const int dimw = IG::dimensionworld;
821 Dune::GeometryType gtface = ig.geometryInInside().type();
822 const int qorder = std::max ( 2 * ( (
int)lfsu.finiteElement().localBasis().order() - 1 ), 0) + superintegration_order;
823 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
827 if( bctype.isDirichlet( ig, rule.begin()->position() ) )
830 const Dune::FieldVector<DF,IG::dimension-1>& face_center =
831 Dune::ReferenceElements<DF,IG::dimension-1>::
832 general(ig.geometry().type()).position(0,0);
834 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(face_center);
837 const Dune::FieldVector<DF,IG::dimension>
838 localcenter = Dune::ReferenceElements<DF,IG::dimension>::general(ig.inside()->type()).position(0,0);
839 typename K::Traits::RangeType tensor(0.0);
840 k.evaluate(*ig.inside(),localcenter,tensor);
843 RF penalty_weight = sigma / pow(ig.geometry().volume(), beta);
846 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
849 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
852 std::vector<JacobianType> js(lfsv.size());
853 lfsv.finiteElement().localBasis().evaluateJacobian(local,js);
856 typename IG::Entity::Geometry::JacobianInverseTransposed jac;
857 jac = ig.inside()->geometry().jacobianInverseTransposed(local);
858 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsv.size());
859 for (
size_t i=0; i<lfsv.size(); i++)
862 jac.umv(js[i][0],gradphi[i]);
866 std::vector<RangeType> phi(lfsv.size());
867 lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
870 Dune::FieldVector<RF,dim> gradu(0.0);
871 for (
size_t i=0; i<lfsu.size(); i++)
873 gradu.axpy(x(lfsu,i),gradphi[i]);
877 std::vector<Dune::FieldVector<RF,dim> > kgradphi(lfsu.size());
878 for (
size_t i=0; i<lfsu.size(); i++)
880 tensor.mv(gradphi[i],kgradphi[i]);
884 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
885 for (
typename LFSU::Traits::SizeType j=0; j<lfsu.size(); j++)
887 for (
typename LFSU::Traits::SizeType i=0; i<lfsu.size(); i++)
890 mat.accumulate( lfsu, i, lfsu, j, (epsilon*(kgradphi[i]*normal)*phi[j] - phi[i]*(kgradphi[j]*normal))*factor );
892 mat.accumulate( lfsu, i, lfsu, j, (penalty_weight*phi[j]*phi[i])*factor );
910 int superintegration_order;
Definition: diffusiondg.hh:59
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusiondg.hh:96
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: diffusiondg.hh:455
Definition: diffusiondg.hh:60
sparsity pattern generator
Definition: pattern.hh:13
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:636
const IG & ig
Definition: constraints.hh:148
static const int dim
Definition: adaptivity.hh:83
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: diffusiondg.hh:573
void jacobian_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: diffusiondg.hh:802
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:157
Implement jacobian_boundary() based on alpha_boundary()
Definition: numericaljacobian.hh:250
Definition: adaptivity.hh:27
Implement alpha_boundary() based on jacobian_boundary()
Definition: numericalresidual.hh:128
sparsity pattern generator
Definition: pattern.hh:29
Default flags for all local operators.
Definition: flags.hh:18
Definition: diffusiondg.hh:54
Definition: diffusiondg.hh:61
Implement jacobian_volume() based on alpha_volume()
Definition: numericaljacobian.hh:31
Definition: diffusiondg.hh:63
Implement alpha_volume() based on jacobian_volume()
Definition: numericalresidual.hh:35
void alpha_boundary(const IG &ig, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusiondg.hh:310
Definition: diffusiondg.hh:58
Definition: diffusiondg.hh:55
Implement jacobian_skeleton() based on alpha_skeleton()
Definition: numericaljacobian.hh:156
Definition: diffusiondg.hh:62
Implement alpha_skeleton() based on jacobian_skeleton()
Definition: numericalresidual.hh:74
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: diffusiondg.hh:413
Definition: diffusiondg.hh:33