2 #ifndef DUNE_PDELAB_DIFFUSION_HH
3 #define DUNE_PDELAB_DIFFUSION_HH
7 #include<dune/common/exceptions.hh>
8 #include<dune/common/fvector.hh>
9 #include<dune/common/static_assert.hh>
10 #include<dune/geometry/type.hh>
11 #include<dune/geometry/referenceelements.hh>
12 #include<dune/geometry/quadraturerules.hh>
38 template<
typename K,
typename A0,
typename F,
typename B,
typename J>
54 Diffusion (
const K& k_,
const A0& a0_,
const F& f_,
const B& bctype_,
const J& j_,
int intorder_=2)
55 : k(k_), a0(a0_), f(f_), bctype(bctype_), j(j_), intorder(intorder_)
59 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
60 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
63 typedef typename LFSU::Traits::FiniteElementType::
64 Traits::LocalBasisType::Traits::DomainFieldType DF;
65 typedef typename LFSU::Traits::FiniteElementType::
66 Traits::LocalBasisType::Traits::RangeFieldType RF;
67 typedef typename LFSU::Traits::FiniteElementType::
68 Traits::LocalBasisType::Traits::JacobianType JacobianType;
69 typedef typename LFSU::Traits::FiniteElementType::
70 Traits::LocalBasisType::Traits::RangeType RangeType;
72 typedef typename LFSU::Traits::SizeType size_type;
75 const int dim = EG::Geometry::dimension;
78 Dune::GeometryType gt = eg.geometry().type();
79 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
82 typename K::Traits::RangeType tensor(0.0);
83 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
84 k.evaluate(eg.entity(),localcenter,tensor);
87 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
90 std::vector<JacobianType> js(lfsu.size());
91 lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
94 const typename EG::Geometry::JacobianInverseTransposed jac =
95 eg.geometry().jacobianInverseTransposed(it->position());
96 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
97 for (size_type i=0; i<lfsu.size(); i++)
100 jac.umv(js[i][0],gradphi[i]);
104 Dune::FieldVector<RF,dim> gradu(0.0);
105 for (size_type i=0; i<lfsu.size(); i++)
106 gradu.axpy(x[i],gradphi[i]);
109 Dune::FieldVector<RF,dim> Kgradu(0.0);
110 tensor.umv(gradu,Kgradu);
113 std::vector<RangeType> phi(lfsu.size());
114 lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
118 for (size_type i=0; i<lfsu.size(); i++)
122 typename A0::Traits::RangeType y;
123 a0.evaluate(eg.entity(),it->position(),y);
126 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
127 for (size_type i=0; i<lfsu.size(); i++)
128 r[i] += ( Kgradu*gradphi[i] + y*u*phi[i] )*factor;
133 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
138 typedef typename LFSU::Traits::FiniteElementType::
139 Traits::LocalBasisType::Traits::DomainFieldType DF;
140 typedef typename LFSU::Traits::FiniteElementType::
141 Traits::LocalBasisType::Traits::RangeFieldType RF;
142 typedef typename LFSU::Traits::FiniteElementType::
143 Traits::LocalBasisType::Traits::JacobianType JacobianType;
144 typedef typename LFSU::Traits::FiniteElementType::
145 Traits::LocalBasisType::Traits::RangeType RangeType;
146 typedef typename LFSU::Traits::SizeType size_type;
149 const int dim = EG::Geometry::dimension;
152 Dune::GeometryType gt = eg.geometry().type();
153 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
156 typename K::Traits::RangeType tensor(0.0);
157 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
158 k.evaluate(eg.entity(),localcenter,tensor);
161 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
164 std::vector<JacobianType> js(lfsu.size());
165 lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
168 const typename EG::Geometry::JacobianInverseTransposed jac =
169 eg.geometry().jacobianInverseTransposed(it->position());
170 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
171 for (size_type i=0; i<lfsu.size(); i++)
174 jac.umv(js[i][0],gradphi[i]);
178 std::vector<Dune::FieldVector<RF,dim> > Kgradphi(lfsu.size());
179 for (size_type i=0; i<lfsu.size(); i++)
180 tensor.mv(gradphi[i],Kgradphi[i]);
183 std::vector<RangeType> phi(lfsu.size());
184 lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
187 typename A0::Traits::RangeType y;
188 a0.evaluate(eg.entity(),it->position(),y);
191 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
192 for (size_type j=0; j<lfsu.size(); j++)
193 for (size_type i=0; i<lfsu.size(); i++)
194 mat(i,j) += ( Kgradphi[j]*gradphi[i] + y*phi[j]*phi[i] )*factor;
199 template<
typename EG,
typename LFSV,
typename R>
203 typedef typename LFSV::Traits::FiniteElementType::
204 Traits::LocalBasisType::Traits::DomainFieldType DF;
205 typedef typename LFSV::Traits::FiniteElementType::
206 Traits::LocalBasisType::Traits::RangeFieldType RF;
207 typedef typename LFSV::Traits::FiniteElementType::
208 Traits::LocalBasisType::Traits::RangeType RangeType;
210 typedef typename LFSV::Traits::SizeType size_type;
213 const int dim = EG::Geometry::dimension;
216 Dune::GeometryType gt = eg.geometry().type();
217 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
220 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
223 std::vector<RangeType> phi(lfsv.size());
224 lfsv.finiteElement().localBasis().evaluateFunction(it->position(),phi);
227 typename F::Traits::RangeType y;
228 f.evaluate(eg.entity(),it->position(),y);
231 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
232 for (size_type i=0; i<lfsv.size(); i++)
233 r[i] -= y*phi[i]*factor;
238 template<
typename IG,
typename LFSV,
typename R>
242 typedef typename LFSV::Traits::FiniteElementType::
243 Traits::LocalBasisType::Traits::DomainFieldType DF;
244 typedef typename LFSV::Traits::FiniteElementType::
245 Traits::LocalBasisType::Traits::RangeFieldType RF;
246 typedef typename LFSV::Traits::FiniteElementType::
247 Traits::LocalBasisType::Traits::RangeType RangeType;
249 typedef typename LFSV::Traits::SizeType size_type;
252 const int dim = IG::dimension;
255 Dune::GeometryType gtface = ig.geometryInInside().type();
256 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
259 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
263 if( bctype.isDirichlet( ig,it->position() ) )
267 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
270 std::vector<RangeType> phi(lfsv.size());
271 lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
274 typename J::Traits::RangeType y;
275 j.evaluate(*(ig.inside()),local,y);
278 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
279 for (size_type i=0; i<lfsv.size(); i++)
280 r[i] += y*phi[i]*factor;
Definition: diffusion.hh:51
Default flags for all local operators.
Definition: flags.hh:18
Definition: diffusion.hh:47
static const int dim
Definition: adaptivity.hh:82
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: diffusion.hh:239
const IG & ig
Definition: common/constraints.hh:146
Definition: diffusion.hh:52
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: diffusion.hh:200
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix< R > &mat) const
Definition: diffusion.hh:134
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusion.hh:60
Definition: diffusion.hh:39
sparsity pattern generator
Definition: pattern.hh:14
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
const EG & eg
Definition: common/constraints.hh:277
Definition: diffusion.hh:50
Diffusion(const K &k_, const A0 &a0_, const F &f_, const B &bctype_, const J &j_, int intorder_=2)
Definition: diffusion.hh:54