2 #ifndef DUNE_PDELAB_CGSTOKES_HH
3 #define DUNE_PDELAB_CGSTOKES_HH
7 #include<dune/common/exceptions.hh>
8 #include<dune/common/fvector.hh>
9 #include<dune/common/static_assert.hh>
11 #include<dune/geometry/type.hh>
12 #include<dune/geometry/quadraturerules.hh>
64 static const bool navier = P::assemble_navier;
80 , _quadrature_order(quadrature_order)
84 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
85 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
88 const int dim = EG::Geometry::dimension;
91 typedef typename LFSU::template Child<0>::Type LFSU_V_PFS;
92 const LFSU_V_PFS& lfsu_v_pfs = lfsu.template child<0>();
94 typedef typename LFSU_V_PFS::template Child<0>::Type LFSU_V;
95 const unsigned int vsize = lfsu_v_pfs.child(0).size();
98 typedef typename LFSU_V::Traits::FiniteElementType::
99 Traits::LocalBasisType::Traits::RangeFieldType RF;
100 typedef typename LFSU_V::Traits::FiniteElementType::
101 Traits::LocalBasisType::Traits::RangeType RT_V;
102 typedef typename LFSU_V::Traits::FiniteElementType::
103 Traits::LocalBasisType::Traits::JacobianType JacobianType_V;
106 typedef typename LFSU::template Child<1>::Type LFSU_P;
107 const LFSU_P& lfsu_p = lfsu.template child<1>();
108 const unsigned int psize = lfsu_p.size();
110 typedef typename LFSU_P::Traits::FiniteElementType::
111 Traits::LocalBasisType::Traits::DomainFieldType DF;
112 typedef typename LFSU_P::Traits::FiniteElementType::
113 Traits::LocalBasisType::Traits::RangeType RT_P;
116 Dune::GeometryType gt = eg.geometry().type();
117 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,_quadrature_order);
120 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it = rule.begin(),
126 std::vector<JacobianType_V> js(vsize);
127 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(it->position(),js);
130 const typename EG::Geometry::JacobianInverseTransposed jac =
131 eg.geometry().jacobianInverseTransposed(it->position());
132 std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
133 for (
size_t i=0; i<vsize; i++)
136 jac.umv(js[i][0],gradphi[i]);
140 std::vector<RT_P> psi(psize);
141 lfsu_p.finiteElement().localBasis().evaluateFunction(it->position(),psi);
144 Dune::FieldVector<RF,dim> vu(0.0);
146 std::vector<RT_V> phi(vsize);
149 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(it->position(),phi);
151 for(
int d=0; d<
dim; ++d)
153 const LFSU_V & lfsu_v = lfsu_v_pfs.child(d);
154 for (
size_t i=0; i<lfsu_v.size(); i++)
155 vu[d] += x(lfsu_v,i) * phi[i];
160 Dune::FieldMatrix<RF,dim,dim> jacu(0.0);
161 for(
int d=0; d<
dim; ++d){
162 const LFSU_V & lfsu_v = lfsu_v_pfs.child(d);
163 for (
size_t i=0; i<lfsu_v.size(); i++)
164 jacu[d].axpy(x(lfsu_v,i),gradphi[i]);
169 for (
size_t i=0; i<lfsu_p.size(); i++)
170 func_p += x(lfsu_p,i) * psi[i];
173 const RF mu = _p.mu(eg,it->position());
174 const RF rho = _p.rho(eg,it->position());
177 const RF factor = it->weight() * eg.geometry().integrationElement(it->position());
179 for(
int d=0; d<
dim; ++d){
181 const LFSU_V & lfsu_v = lfsu_v_pfs.child(d);
184 const RF u_nabla_u = vu * jacu[d];
186 for (
size_t i=0; i<vsize; i++){
189 r.accumulate(lfsu_v,i, mu * (jacu[d] * gradphi[i]) * factor);
193 for(
int dd=0; dd<
dim; ++dd)
194 r.accumulate(lfsu_v,i, mu * (jacu[dd][d] * gradphi[i][dd]) * factor);
197 r.accumulate(lfsu_v,i,- (func_p * gradphi[i][d]) * factor);
201 r.accumulate(lfsu_v,i, rho * u_nabla_u * phi[i] * factor);
208 for (
size_t i=0; i<vsize; i++)
209 for(
int d=0; d<
dim; ++d)
210 divu += x(lfsu_v_pfs.child(d),i) * gradphi[i][d];
213 for (
size_t i=0; i<lfsu_p.size(); i++)
215 r.accumulate(lfsu_p,i, -1.0 * divu * psi[i] * factor);
223 template<
typename EG,
typename LFSV,
typename R>
227 const int dim = EG::Geometry::dimension;
230 typedef typename LFSV::template Child<0>::Type LFSV_V_PFS;
231 const LFSV_V_PFS& lfsv_v_pfs = lfsv.template child<0>();
233 typedef typename LFSV_V_PFS::template Child<0>::Type LFSV_V;
234 const unsigned int vsize = lfsv_v_pfs.child(0).size();
237 typedef typename LFSV_V::Traits::FiniteElementType::
238 Traits::LocalBasisType::Traits::RangeFieldType RF;
239 typedef typename LFSV_V::Traits::FiniteElementType::
240 Traits::LocalBasisType::Traits::RangeType RT_V;
242 typedef typename LFSV::template Child<1>::Type LFSV_P;
243 const LFSV_P& lfsv_p = lfsv.template child<1>();
244 const unsigned int psize = lfsv_p.size();
246 typedef typename LFSV_V::Traits::FiniteElementType::
247 Traits::LocalBasisType::Traits::DomainFieldType DF;
248 typedef typename LFSV_P::Traits::FiniteElementType::
249 Traits::LocalBasisType::Traits::RangeType RT_P;
252 Dune::GeometryType gt = eg.geometry().type();
253 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,_quadrature_order);
256 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it = rule.begin(),
261 std::vector<RT_V> phi(vsize);
262 lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(it->position(),phi);
264 std::vector<RT_P> psi(psize);
265 lfsv_p.finiteElement().localBasis().evaluateFunction(it->position(),psi);
268 const Dune::FieldVector<RF,dim> f1 = _p.f(eg,it->position());
271 const RF factor = it->weight() * eg.geometry().integrationElement(it->position());
273 for(
int d=0; d<
dim; ++d){
275 const LFSV_V & lfsv_v = lfsv_v_pfs.child(d);
277 for (
size_t i=0; i<vsize; i++)
280 r.accumulate(lfsv_v,i, -f1[d]*phi[i] * factor);
285 const RF g2 = _p.g2(eg,it->position());
288 for (
size_t i=0; i<lfsv_p.size(); i++)
290 r.accumulate(lfsv_p,i, g2 * psi[i] * factor);
298 template<
typename IG,
typename LFSV,
typename R>
302 static const unsigned int dim = IG::Geometry::dimension;
303 static const unsigned int dimw = IG::Geometry::dimensionworld;
306 typedef typename LFSV::template Child<0>::Type LFSV_V_PFS;
307 const LFSV_V_PFS& lfsv_v_pfs = lfsv.template child<0>();
309 typedef typename LFSV_V_PFS::template Child<0>::Type LFSV_V;
310 const unsigned int vsize = lfsv_v_pfs.child(0).size();
312 typedef typename LFSV_V::Traits::FiniteElementType::
313 Traits::LocalBasisType::Traits::RangeType RT_V;
316 typedef typename LFSV_V::Traits::FiniteElementType::
317 Traits::LocalBasisType::Traits::RangeFieldType RF;
320 typedef typename LFSV_V::Traits::FiniteElementType::
321 Traits::LocalBasisType::Traits::DomainFieldType DF;
324 Dune::GeometryType gtface = ig.geometryInInside().type();
325 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,_quadrature_order);
328 for (
typename Dune::QuadratureRule<DF,dim-1>::const_iterator it = rule.begin(),
334 typename BC::Type bctype = _p.bctype(ig,it->position());
341 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(it->position());
344 std::vector<RT_V> phi(vsize);
345 lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(local,phi);
347 const RF factor = it->weight() * ig.geometry().integrationElement(it->position());
348 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(it->position());
351 const Dune::FieldVector<DF,dimw> neumann_stress = _p.j(ig,it->position(),normal);
353 for(
unsigned int d=0; d<
dim; ++d)
356 const LFSV_V & lfsv_v = lfsv_v_pfs.child(d);
358 for (
size_t i=0; i<vsize; i++)
360 r.accumulate(lfsv_v,i, neumann_stress[d] * phi[i] * factor);
368 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
373 const int dim = EG::Geometry::dimension;
377 typedef typename LFSU::template Child<0>::Type LFSU_V_PFS;
378 const LFSU_V_PFS& lfsu_v_pfs = lfsu.template child<0>();
379 const unsigned int vsize = lfsu_v_pfs.child(0).size();
381 typedef typename LFSU_V_PFS::template Child<0>::Type LFSU_V;
384 typedef typename LFSU_V::Traits::FiniteElementType::
385 Traits::LocalBasisType::Traits::RangeFieldType RF;
386 typedef typename LFSU_V::Traits::FiniteElementType::
387 Traits::LocalBasisType::Traits::RangeType RT_V;
388 typedef typename LFSU_V::Traits::FiniteElementType::
389 Traits::LocalBasisType::Traits::JacobianType JacobianType_V;
392 typedef typename LFSU::template Child<1>::Type LFSU_P;
393 const LFSU_P& lfsu_p = lfsu.template child<1>();
394 const unsigned int psize = lfsu_p.size();
396 typedef typename LFSU_P::Traits::FiniteElementType::
397 Traits::LocalBasisType::Traits::DomainFieldType DF;
399 typedef typename LFSU_P::Traits::FiniteElementType::
400 Traits::LocalBasisType::Traits::RangeType RT_P;
403 Dune::GeometryType gt = eg.geometry().type();
404 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,_quadrature_order);
407 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it = rule.begin(),
413 std::vector<JacobianType_V> js(vsize);
414 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(it->position(),js);
417 const typename EG::Geometry::JacobianInverseTransposed jac =
418 eg.geometry().jacobianInverseTransposed(it->position());
419 std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
420 for (
size_t i=0; i<vsize; i++)
423 jac.umv(js[i][0],gradphi[i]);
427 std::vector<RT_P> psi(psize);
428 lfsu_p.finiteElement().localBasis().evaluateFunction(it->position(),psi);
431 std::vector<RT_V> phi(vsize);
432 Dune::FieldVector<RF,dim> vu(0.0);
434 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(it->position(),phi);
436 for(
int d = 0; d <
dim; ++d){
437 const LFSU_V & lfsv_v = lfsu_v_pfs.child(d);
438 for(
size_t l = 0; l < vsize; ++l)
439 vu[d] += x(lfsv_v,l) * phi[l];
444 const RF mu = _p.mu(eg,it->position());
445 const RF rho = _p.rho(eg,it->position());
447 const RF factor = it->weight() * eg.geometry().integrationElement(it->position());
449 for(
int d=0; d<
dim; ++d){
451 const LFSU_V & lfsv_v = lfsu_v_pfs.child(d);
452 const LFSU_V & lfsu_v = lfsv_v;
455 Dune::FieldVector<RF,dim> gradu_d(0.0);
457 for(
size_t l =0; l < vsize; ++l)
458 gradu_d.axpy(x(lfsv_v,l), gradphi[l]);
460 for (
size_t i=0; i<lfsv_v.size(); i++){
463 for (
size_t j=0; j<lfsv_v.size(); j++){
464 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[i] * gradphi[j]) * factor);
468 for(
int dd=0; dd<
dim; ++dd){
469 const LFSU_V & lfsu_v = lfsu_v_pfs.child(dd);
470 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[j][d] * gradphi[i][dd]) * factor);
476 for (
size_t j=0; j<lfsu_p.size(); j++)
477 mat.accumulate(lfsv_v,i,lfsu_p,j, - (gradphi[i][d] * psi[j]) * factor);
480 for(
int k =0; k <
dim; ++k){
481 const LFSU_V & lfsu_v = lfsu_v_pfs.child(k);
483 const RF pre_factor = factor * rho * gradu_d[k] * phi[i];
485 for(
size_t j=0; j< lfsu_v.size(); ++j)
486 mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * phi[j]);
491 const RF pre_factor = factor * rho * phi[i];
492 for(
size_t j=0; j< lfsu_v.size(); ++j)
493 mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * (vu * gradphi[j]));
498 for (
size_t i=0; i<lfsu_p.size(); i++){
499 for (
size_t j=0; j<lfsu_v.size(); j++)
500 mat.accumulate(lfsu_p,i,lfsu_v,j, - (gradphi[j][d] * psi[i]) * factor);
510 const std::size_t _quadrature_order;
521 template<
typename P >
536 : p(p_), intorder(intorder_)
540 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
541 void alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
543 typedef typename LFSV::template Child<0>::Type LFSV_PFS_V;
544 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<0>();
546 for(
unsigned int i=0; i<LFSV_PFS_V::CHILDREN; ++i)
548 scalar_alpha_volume(eg,lfsv_pfs_v.child(i),x,lfsv_pfs_v.child(i),r);
553 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
557 typedef typename LFSV::template Child<0>::Type LFSV_PFS_V;
558 const LFSV_PFS_V& lfsv_pfs_v = lfsv.template child<0>();
560 for(
unsigned int i=0; i<LFSV_PFS_V::CHILDREN; ++i)
562 scalar_jacobian_volume(eg,lfsv_pfs_v.child(i),x,lfsv_pfs_v.child(i),mat);
567 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
568 void scalar_alpha_volume (
const EG&
eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
573 typedef FiniteElementInterfaceSwitch<
574 typename LFSU::Traits::FiniteElementType
576 typedef BasisInterfaceSwitch<
577 typename FESwitch::Basis
581 typedef typename BasisSwitch::DomainField DF;
582 typedef typename BasisSwitch::RangeField RF;
583 typedef typename BasisSwitch::Range RangeType;
585 typedef typename LFSU::Traits::SizeType size_type;
588 const int dim = EG::Geometry::dimension;
591 Dune::GeometryType gt = eg.geometry().type();
592 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
595 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
598 std::vector<RangeType> phi(lfsu.size());
599 FESwitch::basis(lfsu.finiteElement()).evaluateFunction(it->position(),phi);
601 RF rho = p.rho(eg,it->position());
604 for (size_type i=0; i<lfsu.size(); i++)
605 u += x(lfsu,i)*phi[i];
608 RF factor = it->weight() * rho * eg.geometry().integrationElement(it->position());
610 for (size_type i=0; i<lfsu.size(); i++)
611 r.accumulate(lfsv,i, u*phi[i]*factor);
615 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
616 void scalar_jacobian_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
621 typedef FiniteElementInterfaceSwitch<
622 typename LFSU::Traits::FiniteElementType
624 typedef BasisInterfaceSwitch<
625 typename FESwitch::Basis
629 typedef typename BasisSwitch::DomainField DF;
630 typedef typename BasisSwitch::RangeField RF;
631 typedef typename BasisSwitch::Range RangeType;
632 typedef typename LFSU::Traits::SizeType size_type;
635 const int dim = EG::Geometry::dimension;
638 Dune::GeometryType gt = eg.geometry().type();
639 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,intorder);
642 for (
typename Dune::QuadratureRule<DF,dim>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
645 std::vector<RangeType> phi(lfsu.size());
646 FESwitch::basis(lfsu.finiteElement()).evaluateFunction(it->position(),phi);
649 RF rho = p.rho(eg,it->position());
650 RF factor = it->weight() * rho * eg.geometry().integrationElement(it->position());
651 for (size_type j=0; j<lfsu.size(); j++)
652 for (size_type i=0; i<lfsu.size(); i++)
653 mat.accumulate(lfsv,i,lfsu,j, phi[j]*phi[i]*factor);
Type
Definition: stokesparameter.hh:33
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: cg_stokes.hh:85
Definition: cg_stokes.hh:68
Default flags for all local operators.
Definition: flags.hh:18
A local operator for the Navier-Stokes equations.
Definition: cg_stokes.hh:55
NavierStokesMass(const P &p_, int intorder_=4)
Definition: cg_stokes.hh:535
static const bool navier
Definition: cg_stokes.hh:64
P PhysicalParameters
Definition: cg_stokes.hh:75
Definition: cg_stokes.hh:533
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: cg_stokes.hh:541
TaylorHoodNavierStokes(const PhysicalParameters &p, std::size_t quadrature_order=4)
Definition: cg_stokes.hh:77
static const int dim
Definition: adaptivity.hh:82
StokesBoundaryCondition BC
Boundary condition indicator type.
Definition: cg_stokes.hh:62
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: cg_stokes.hh:299
Definition: stokesparameter.hh:32
Definition: cg_stokes.hh:72
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: cg_stokes.hh:224
Definition: cg_stokes.hh:530
const IG & ig
Definition: common/constraints.hh:146
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: cg_stokes.hh:554
Definition: cg_stokes.hh:522
Definition: cg_stokes.hh:73
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
sparsity pattern generator
Definition: pattern.hh:14
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
static const bool full_tensor
Definition: cg_stokes.hh:65
Definition: stokesparameter.hh:36
const EG & eg
Definition: common/constraints.hh:277
Definition: cg_stokes.hh:71
R p(int i, D x)
Definition: qkdg.hh:62
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: cg_stokes.hh:369