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/geometry/type.hh>
10 #include<dune/geometry/referenceelements.hh>
11 #include<dune/geometry/quadraturerules.hh>
37 template<
typename K,
typename A0,
typename F,
typename B,
typename J>
53 Diffusion (
const K& k_,
const A0& a0_,
const F& f_,
const B& bctype_,
const J& j_,
int intorder_=2)
54 : k(k_), a0(a0_), f(f_), bctype(bctype_), j(j_), intorder(intorder_)
58 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
59 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
62 typedef typename LFSU::Traits::FiniteElementType::
63 Traits::LocalBasisType::Traits::DomainFieldType DF;
64 typedef typename LFSU::Traits::FiniteElementType::
65 Traits::LocalBasisType::Traits::RangeFieldType RF;
66 typedef typename LFSU::Traits::FiniteElementType::
67 Traits::LocalBasisType::Traits::JacobianType JacobianType;
68 typedef typename LFSU::Traits::FiniteElementType::
69 Traits::LocalBasisType::Traits::RangeType RangeType;
71 typedef typename LFSU::Traits::SizeType size_type;
74 const int dim = EG::Geometry::dimension;
77 Dune::GeometryType gt = eg.geometry().type();
78 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
81 typename K::Traits::RangeType tensor(0.0);
82 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
83 k.evaluate(eg.entity(),localcenter,tensor);
86 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
89 std::vector<JacobianType> js(lfsu.size());
90 lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
93 const typename EG::Geometry::JacobianInverseTransposed jac =
94 eg.geometry().jacobianInverseTransposed(it->position());
95 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
96 for (size_type i=0; i<lfsu.size(); i++)
99 jac.umv(js[i][0],gradphi[i]);
103 Dune::FieldVector<RF,dim> gradu(0.0);
104 for (size_type i=0; i<lfsu.size(); i++)
105 gradu.axpy(x[i],gradphi[i]);
108 Dune::FieldVector<RF,dim> Kgradu(0.0);
109 tensor.umv(gradu,Kgradu);
112 std::vector<RangeType> phi(lfsu.size());
113 lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
117 for (size_type i=0; i<lfsu.size(); i++)
121 typename A0::Traits::RangeType y;
122 a0.evaluate(eg.entity(),it->position(),y);
125 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
126 for (size_type i=0; i<lfsu.size(); i++)
127 r[i] += ( Kgradu*gradphi[i] + y*u*phi[i] )*factor;
132 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
137 typedef typename LFSU::Traits::FiniteElementType::
138 Traits::LocalBasisType::Traits::DomainFieldType DF;
139 typedef typename LFSU::Traits::FiniteElementType::
140 Traits::LocalBasisType::Traits::RangeFieldType RF;
141 typedef typename LFSU::Traits::FiniteElementType::
142 Traits::LocalBasisType::Traits::JacobianType JacobianType;
143 typedef typename LFSU::Traits::FiniteElementType::
144 Traits::LocalBasisType::Traits::RangeType RangeType;
145 typedef typename LFSU::Traits::SizeType size_type;
148 const int dim = EG::Geometry::dimension;
151 Dune::GeometryType gt = eg.geometry().type();
152 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
155 typename K::Traits::RangeType tensor(0.0);
156 Dune::FieldVector<DF,dim> localcenter = Dune::ReferenceElements<DF,dim>::general(gt).position(0,0);
157 k.evaluate(eg.entity(),localcenter,tensor);
160 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
163 std::vector<JacobianType> js(lfsu.size());
164 lfsu.finiteElement().localBasis().evaluateJacobian(it->position(),js);
167 const typename EG::Geometry::JacobianInverseTransposed jac =
168 eg.geometry().jacobianInverseTransposed(it->position());
169 std::vector<Dune::FieldVector<RF,dim> > gradphi(lfsu.size());
170 for (size_type i=0; i<lfsu.size(); i++)
173 jac.umv(js[i][0],gradphi[i]);
177 std::vector<Dune::FieldVector<RF,dim> > Kgradphi(lfsu.size());
178 for (size_type i=0; i<lfsu.size(); i++)
179 tensor.mv(gradphi[i],Kgradphi[i]);
182 std::vector<RangeType> phi(lfsu.size());
183 lfsu.finiteElement().localBasis().evaluateFunction(it->position(),phi);
186 typename A0::Traits::RangeType y;
187 a0.evaluate(eg.entity(),it->position(),y);
190 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
191 for (size_type j=0; j<lfsu.size(); j++)
192 for (size_type i=0; i<lfsu.size(); i++)
193 mat(i,j) += ( Kgradphi[j]*gradphi[i] + y*phi[j]*phi[i] )*factor;
198 template<
typename EG,
typename LFSV,
typename R>
202 typedef typename LFSV::Traits::FiniteElementType::
203 Traits::LocalBasisType::Traits::DomainFieldType DF;
204 typedef typename LFSV::Traits::FiniteElementType::
205 Traits::LocalBasisType::Traits::RangeFieldType RF;
206 typedef typename LFSV::Traits::FiniteElementType::
207 Traits::LocalBasisType::Traits::RangeType RangeType;
209 typedef typename LFSV::Traits::SizeType size_type;
212 const int dim = EG::Geometry::dimension;
215 Dune::GeometryType gt = eg.geometry().type();
216 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
219 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
222 std::vector<RangeType> phi(lfsv.size());
223 lfsv.finiteElement().localBasis().evaluateFunction(it->position(),phi);
226 typename F::Traits::RangeType y;
227 f.evaluate(eg.entity(),it->position(),y);
230 RF factor = it->weight() * eg.geometry().integrationElement(it->position());
231 for (size_type i=0; i<lfsv.size(); i++)
232 r[i] -= y*phi[i]*factor;
237 template<
typename IG,
typename LFSV,
typename R>
241 typedef typename LFSV::Traits::FiniteElementType::
242 Traits::LocalBasisType::Traits::DomainFieldType DF;
243 typedef typename LFSV::Traits::FiniteElementType::
244 Traits::LocalBasisType::Traits::RangeFieldType RF;
245 typedef typename LFSV::Traits::FiniteElementType::
246 Traits::LocalBasisType::Traits::RangeType RangeType;
248 typedef typename LFSV::Traits::SizeType size_type;
251 const int dim = IG::dimension;
254 Dune::GeometryType gtface = ig.geometryInInside().type();
255 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,intorder);
258 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
262 if( bctype.isDirichlet( ig,it->position() ) )
266 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
269 std::vector<RangeType> phi(lfsv.size());
270 lfsv.finiteElement().localBasis().evaluateFunction(local,phi);
273 typename J::Traits::RangeType y;
274 j.evaluate(*(ig.inside()),local,y);
277 RF factor = it->weight()*ig.geometry().integrationElement(it->position());
278 for (size_type i=0; i<lfsv.size(); i++)
279 r[i] += y*phi[i]*factor;
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: diffusion.hh:199
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
const IG & ig
Definition: constraints.hh:147
Definition: diffusion.hh:46
Definition: diffusion.hh:49
Default flags for all local operators.
Definition: flags.hh:18
Definition: diffusion.hh:38
sparsity pattern generator
Definition: pattern.hh:13
Diffusion(const K &k_, const A0 &a0_, const F &f_, const B &bctype_, const J &j_, int intorder_=2)
Definition: diffusion.hh:53
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: diffusion.hh:59
Definition: diffusion.hh:50
Definition: diffusion.hh:51
Definition: adaptivity.hh:27
static const int dim
Definition: adaptivity.hh:83
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, LocalMatrix< R > &mat) const
Definition: diffusion.hh:133
const EG & eg
Definition: constraints.hh:280
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: diffusion.hh:238