2 #ifndef DUNE_PDELAB_LOCALOPERATOR_TAYLORHOODNAVIERSTOKES_HH
3 #define DUNE_PDELAB_LOCALOPERATOR_TAYLORHOODNAVIERSTOKES_HH
7 #include<dune/common/exceptions.hh>
8 #include<dune/common/fvector.hh>
10 #include<dune/geometry/type.hh>
11 #include<dune/geometry/quadraturerules.hh>
64 static const bool navier = P::assemble_navier;
80 , superintegration_order(superintegration_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 int v_order = lfsu_v_pfs.child(0).finiteElement().localBasis().order();
118 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
119 const int jac_order = gt.isSimplex() ? 0 : 1;
120 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
121 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
124 for (
const auto& ip : rule)
127 std::vector<JacobianType_V> js(vsize);
128 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(ip.position(),js);
131 const typename EG::Geometry::JacobianInverseTransposed jac =
132 eg.geometry().jacobianInverseTransposed(ip.position());
133 std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
134 for (
size_t i=0; i<vsize; i++)
137 jac.umv(js[i][0],gradphi[i]);
141 std::vector<RT_P> psi(psize);
142 lfsu_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
145 Dune::FieldVector<RF,dim> vu(0.0);
147 std::vector<RT_V> phi(vsize);
150 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
152 for(
int d=0; d<
dim; ++d)
154 const LFSU_V & lfsu_v = lfsu_v_pfs.child(d);
155 for (
size_t i=0; i<lfsu_v.size(); i++)
156 vu[d] += x(lfsu_v,i) * phi[i];
161 Dune::FieldMatrix<RF,dim,dim> jacu(0.0);
162 for(
int d=0; d<
dim; ++d){
163 const LFSU_V & lfsu_v = lfsu_v_pfs.child(d);
164 for (
size_t i=0; i<lfsu_v.size(); i++)
165 jacu[d].axpy(x(lfsu_v,i),gradphi[i]);
170 for (
size_t i=0; i<lfsu_p.size(); i++)
171 func_p += x(lfsu_p,i) * psi[i];
174 const RF mu = _p.mu(eg,ip.position());
175 const RF rho = _p.rho(eg,ip.position());
178 const RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
180 for(
int d=0; d<
dim; ++d){
182 const LFSU_V & lfsu_v = lfsu_v_pfs.child(d);
185 const RF u_nabla_u = vu * jacu[d];
187 for (
size_t i=0; i<vsize; i++){
190 r.accumulate(lfsu_v,i, mu * (jacu[d] * gradphi[i]) * factor);
194 for(
int dd=0; dd<
dim; ++dd)
195 r.accumulate(lfsu_v,i, mu * (jacu[dd][d] * gradphi[i][dd]) * factor);
198 r.accumulate(lfsu_v,i,- (func_p * gradphi[i][d]) * factor);
202 r.accumulate(lfsu_v,i, rho * u_nabla_u * phi[i] * factor);
208 for (
size_t i=0; i<psize; i++)
209 for(
int d=0; d<
dim; ++d)
211 r.accumulate(lfsu_p,i, -1.0 * jacu[d][d] * psi[i] * factor);
218 template<
typename EG,
typename LFSV,
typename R>
222 const int dim = EG::Geometry::dimension;
225 typedef typename LFSV::template Child<0>::Type LFSV_V_PFS;
226 const LFSV_V_PFS& lfsv_v_pfs = lfsv.template child<0>();
228 typedef typename LFSV_V_PFS::template Child<0>::Type LFSV_V;
229 const unsigned int vsize = lfsv_v_pfs.child(0).size();
232 typedef typename LFSV_V::Traits::FiniteElementType::
233 Traits::LocalBasisType::Traits::RangeFieldType RF;
234 typedef typename LFSV_V::Traits::FiniteElementType::
235 Traits::LocalBasisType::Traits::RangeType RT_V;
237 typedef typename LFSV::template Child<1>::Type LFSV_P;
238 const LFSV_P& lfsv_p = lfsv.template child<1>();
239 const unsigned int psize = lfsv_p.size();
241 typedef typename LFSV_V::Traits::FiniteElementType::
242 Traits::LocalBasisType::Traits::DomainFieldType DF;
243 typedef typename LFSV_P::Traits::FiniteElementType::
244 Traits::LocalBasisType::Traits::RangeType RT_P;
247 Dune::GeometryType gt = eg.geometry().type();
248 const int v_order = lfsv_v_pfs.child(0).finiteElement().localBasis().order();
249 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
250 const int qorder = 2*v_order + det_jac_order + superintegration_order;
251 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
254 for (
const auto& ip : rule)
256 std::vector<RT_V> phi(vsize);
257 lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
259 std::vector<RT_P> psi(psize);
260 lfsv_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
263 const Dune::FieldVector<RF,dim> f1 = _p.f(eg,ip.position());
266 const RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
268 for(
int d=0; d<
dim; ++d){
270 const LFSV_V & lfsv_v = lfsv_v_pfs.child(d);
272 for (
size_t i=0; i<vsize; i++)
275 r.accumulate(lfsv_v,i, -f1[d]*phi[i] * factor);
280 const RF g2 = _p.g2(eg,ip.position());
283 for (
size_t i=0; i<psize; i++)
285 r.accumulate(lfsv_p,i, g2 * psi[i] * factor);
293 template<
typename IG,
typename LFSV,
typename R>
297 static const unsigned int dim = IG::Geometry::dimension;
298 static const unsigned int dimw = IG::Geometry::dimensionworld;
301 typedef typename LFSV::template Child<0>::Type LFSV_V_PFS;
302 const LFSV_V_PFS& lfsv_v_pfs = lfsv.template child<0>();
304 typedef typename LFSV_V_PFS::template Child<0>::Type LFSV_V;
305 const unsigned int vsize = lfsv_v_pfs.child(0).size();
307 typedef typename LFSV_V::Traits::FiniteElementType::
308 Traits::LocalBasisType::Traits::RangeType RT_V;
311 typedef typename LFSV_V::Traits::FiniteElementType::
312 Traits::LocalBasisType::Traits::RangeFieldType RF;
315 typedef typename LFSV_V::Traits::FiniteElementType::
316 Traits::LocalBasisType::Traits::DomainFieldType DF;
319 Dune::GeometryType gtface = ig.geometryInInside().type();
320 const int v_order = lfsv_v_pfs.child(0).finiteElement().localBasis().order();
321 const int det_jac_order = gtface.isSimplex() ? 0 : (dim-2);
322 const int jac_order = gtface.isSimplex() ? 0 : 1;
323 const int qorder = 2*v_order + det_jac_order + jac_order + superintegration_order;
324 const Dune::QuadratureRule<DF,dim-1>& rule = Dune::QuadratureRules<DF,dim-1>::rule(gtface,qorder);
327 for (
const auto& ip : rule)
330 typename BC::Type bctype = _p.bctype(ig,ip.position());
337 Dune::FieldVector<DF,dim> local = ig.geometryInInside().global(ip.position());
340 std::vector<RT_V> phi(vsize);
341 lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(local,phi);
343 const RF factor = ip.weight() * ig.geometry().integrationElement(ip.position());
344 const Dune::FieldVector<DF,dimw> normal = ig.unitOuterNormal(ip.position());
347 const Dune::FieldVector<DF,dimw> neumann_stress = _p.j(ig,ip.position(),normal);
349 for(
unsigned int d=0; d<
dim; ++d)
352 const LFSV_V & lfsv_v = lfsv_v_pfs.child(d);
354 for (
size_t i=0; i<vsize; i++)
356 r.accumulate(lfsv_v,i, neumann_stress[d] * phi[i] * factor);
364 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
369 const int dim = EG::Geometry::dimension;
373 typedef typename LFSU::template Child<0>::Type LFSU_V_PFS;
374 const LFSU_V_PFS& lfsu_v_pfs = lfsu.template child<0>();
375 const unsigned int vsize = lfsu_v_pfs.child(0).size();
377 typedef typename LFSU_V_PFS::template Child<0>::Type LFSU_V;
380 typedef typename LFSU_V::Traits::FiniteElementType::
381 Traits::LocalBasisType::Traits::RangeFieldType RF;
382 typedef typename LFSU_V::Traits::FiniteElementType::
383 Traits::LocalBasisType::Traits::RangeType RT_V;
384 typedef typename LFSU_V::Traits::FiniteElementType::
385 Traits::LocalBasisType::Traits::JacobianType JacobianType_V;
388 typedef typename LFSU::template Child<1>::Type LFSU_P;
389 const LFSU_P& lfsu_p = lfsu.template child<1>();
390 const unsigned int psize = lfsu_p.size();
392 typedef typename LFSU_P::Traits::FiniteElementType::
393 Traits::LocalBasisType::Traits::DomainFieldType DF;
395 typedef typename LFSU_P::Traits::FiniteElementType::
396 Traits::LocalBasisType::Traits::RangeType RT_P;
399 Dune::GeometryType gt = eg.geometry().type();
400 const int v_order = lfsu_v_pfs.child(0).finiteElement().localBasis().order();
401 const int det_jac_order = gt.isSimplex() ? 0 : (dim-1);
402 const int jac_order = gt.isSimplex() ? 0 : 1;
403 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
404 const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(gt,qorder);
407 for (
const auto& ip : rule)
410 std::vector<JacobianType_V> js(vsize);
411 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(ip.position(),js);
414 const typename EG::Geometry::JacobianInverseTransposed jac =
415 eg.geometry().jacobianInverseTransposed(ip.position());
416 std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
417 for (
size_t i=0; i<vsize; i++)
420 jac.umv(js[i][0],gradphi[i]);
424 std::vector<RT_P> psi(psize);
425 lfsu_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
428 std::vector<RT_V> phi(vsize);
429 Dune::FieldVector<RF,dim> vu(0.0);
431 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
433 for(
int d = 0; d <
dim; ++d){
434 const LFSU_V & lfsv_v = lfsu_v_pfs.child(d);
435 for(
size_t l = 0; l < vsize; ++l)
436 vu[d] += x(lfsv_v,l) * phi[l];
441 const RF mu = _p.mu(eg,ip.position());
442 const RF rho = _p.rho(eg,ip.position());
444 const RF factor = ip.weight() * eg.geometry().integrationElement(ip.position());
446 for(
int d=0; d<
dim; ++d){
448 const LFSU_V & lfsv_v = lfsu_v_pfs.child(d);
449 const LFSU_V & lfsu_v = lfsv_v;
452 Dune::FieldVector<RF,dim> gradu_d(0.0);
454 for(
size_t l =0; l < vsize; ++l)
455 gradu_d.axpy(x(lfsv_v,l), gradphi[l]);
457 for (
size_t i=0; i<lfsv_v.size(); i++){
460 for (
size_t j=0; j<lfsv_v.size(); j++){
461 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[i] * gradphi[j]) * factor);
465 for(
int dd=0; dd<
dim; ++dd){
466 const LFSU_V & lfsu_v = lfsu_v_pfs.child(dd);
467 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[j][d] * gradphi[i][dd]) * factor);
473 for (
size_t j=0; j<psize; j++)
474 mat.accumulate(lfsv_v,i,lfsu_p,j, - (gradphi[i][d] * psi[j]) * factor);
477 for(
int k =0; k <
dim; ++k){
478 const LFSU_V & lfsu_v = lfsu_v_pfs.child(k);
480 const RF pre_factor = factor * rho * gradu_d[k] * phi[i];
482 for(
size_t j=0; j< lfsu_v.size(); ++j)
483 mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * phi[j]);
488 const RF pre_factor = factor * rho * phi[i];
489 for(
size_t j=0; j< lfsu_v.size(); ++j)
490 mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * (vu * gradphi[j]));
495 for (
size_t i=0; i<psize; i++){
496 for (
size_t j=0; j<lfsu_v.size(); j++)
497 mat.accumulate(lfsu_p,i,lfsu_v,j, - (gradphi[j][d] * psi[i]) * factor);
507 const int superintegration_order;
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: taylorhoodnavierstokes.hh:365
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:85
TaylorHoodNavierStokes(const PhysicalParameters &p, int superintegration_order_=0)
Definition: taylorhoodnavierstokes.hh:77
static const bool navier
Definition: taylorhoodnavierstokes.hh:64
Definition: taylorhoodnavierstokes.hh:73
Definition: stokesparameter.hh:36
P PhysicalParameters
Definition: taylorhoodnavierstokes.hh:75
StokesBoundaryCondition BC
Boundary condition indicator type.
Definition: taylorhoodnavierstokes.hh:62
Definition: taylorhoodnavierstokes.hh:71
const IG & ig
Definition: constraints.hh:147
Default flags for all local operators.
Definition: flags.hh:18
sparsity pattern generator
Definition: pattern.hh:13
Definition: stokesparameter.hh:32
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:294
Definition: taylorhoodnavierstokes.hh:68
Definition: adaptivity.hh:27
Type
Definition: stokesparameter.hh:33
static const int dim
Definition: adaptivity.hh:83
A local operator for the Navier-Stokes equations.
Definition: taylorhoodnavierstokes.hh:55
const P & p
Definition: constraints.hh:146
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:219
static const bool full_tensor
Definition: taylorhoodnavierstokes.hh:65
Definition: taylorhoodnavierstokes.hh:72
const EG & eg
Definition: constraints.hh:280
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89