2#ifndef DUNE_PDELAB_LOCALOPERATOR_TAYLORHOODNAVIERSTOKES_HH
3#define DUNE_PDELAB_LOCALOPERATOR_TAYLORHOODNAVIERSTOKES_HH
11#include<dune/geometry/referenceelements.hh>
13#include<dune/pdelab/common/quadraturerules.hh>
16#include"defaultimp.hh"
21#include"stokesparameter.hh"
22#include"navierstokesmass.hh"
62 using Real =
typename InstatBase::RealType;
64 static const bool navier = P::assemble_navier;
65 static const bool full_tensor = P::assemble_full_tensor;
68 enum { doPatternVolume =
true };
71 enum { doAlphaVolume =
true };
72 enum { doLambdaVolume =
true };
73 enum { doLambdaBoundary =
true };
75 using PhysicalParameters = P;
77 TaylorHoodNavierStokes (PhysicalParameters& p,
int superintegration_order_ = 0)
80 , superintegration_order(superintegration_order_)
91 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
92 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const
96 using LFSU_V_PFS = TypeTree::Child<LFSU,_0>;
97 using LFSU_V = TypeTree::Child<LFSU_V_PFS,_0>;
98 using LFSU_P = TypeTree::Child<LFSU,_1>;
99 using RF =
typename LFSU_V::Traits::FiniteElementType::
100 Traits::LocalBasisType::Traits::RangeFieldType;
101 using RT_V =
typename LFSU_V::Traits::FiniteElementType::
102 Traits::LocalBasisType::Traits::RangeType;
103 using JacobianType_V =
typename LFSU_V::Traits::FiniteElementType::
104 Traits::LocalBasisType::Traits::JacobianType;
105 using RT_P =
typename LFSU_P::Traits::FiniteElementType::
106 Traits::LocalBasisType::Traits::RangeType;
109 const auto& lfsu_v_pfs =
child(lfsu,
_0);
110 const unsigned int vsize = lfsu_v_pfs.child(0).size();
111 const auto& lfsu_p =
child(lfsu,
_1);
112 const unsigned int psize = lfsu_p.size();
115 const int dim = EG::Geometry::mydimension;
118 auto geo = eg.geometry();
121 const int v_order = lfsu_v_pfs.child(0).finiteElement().localBasis().order();
122 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
123 const int jac_order = geo.type().isSimplex() ? 0 : 1;
124 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
127 typename EG::Geometry::JacobianInverseTransposed jac;
128 std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
129 std::vector<RT_P> psi(psize);
131 std::vector<RT_V> phi(vsize);
135 for (
const auto& ip : quadratureRule(geo,qorder))
138 std::vector<JacobianType_V> js(vsize);
139 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(ip.position(),js);
142 jac = geo.jacobianInverseTransposed(ip.position());
143 for (
size_t i=0; i<vsize; i++)
146 jac.umv(js[i][0],gradphi[i]);
150 lfsu_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
155 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
157 for(
int d=0; d<dim; ++d)
160 const auto& lfsu_v = lfsu_v_pfs.child(d);
161 for (
size_t i=0; i<lfsu_v.size(); i++)
162 vu[d] += x(lfsu_v,i) * phi[i];
167 for(
int d=0; d<dim; ++d){
169 const auto& lfsu_v = lfsu_v_pfs.child(d);
170 for (
size_t i=0; i<lfsu_v.size(); i++)
171 jacu[d].axpy(x(lfsu_v,i),gradphi[i]);
176 for (
size_t i=0; i<lfsu_p.size(); i++)
177 func_p += x(lfsu_p,i) * psi[i];
180 const auto mu = _p.mu(eg,ip.position());
181 const auto rho = _p.rho(eg,ip.position());
184 const auto factor = ip.weight() * geo.integrationElement(ip.position());
186 for(
int d=0; d<dim; ++d){
188 const auto& lfsu_v = lfsu_v_pfs.child(d);
191 const auto u_nabla_u = vu * jacu[d];
193 for (
size_t i=0; i<vsize; i++){
196 r.accumulate(lfsu_v,i, mu * (jacu[d] * gradphi[i]) * factor);
200 for(
int dd=0; dd<dim; ++dd)
201 r.accumulate(lfsu_v,i, mu * (jacu[dd][d] * gradphi[i][dd]) * factor);
204 r.accumulate(lfsu_v,i,- (func_p * gradphi[i][d]) * factor);
208 r.accumulate(lfsu_v,i, rho * u_nabla_u * phi[i] * factor);
214 for (
size_t i=0; i<psize; i++)
215 for(
int d=0; d<dim; ++d)
217 r.accumulate(lfsu_p,i, -1.0 * jacu[d][d] * psi[i] * factor);
224 template<
typename EG,
typename LFSV,
typename R>
225 void lambda_volume (
const EG& eg,
const LFSV& lfsv, R& r)
const
229 using LFSV_V_PFS = TypeTree::Child<LFSV,_0>;
230 using LFSV_V = TypeTree::Child<LFSV_V_PFS,_0>;
231 using LFSV_P = TypeTree::Child<LFSV,_1>;
232 using RT_V =
typename LFSV_V::Traits::FiniteElementType::
233 Traits::LocalBasisType::Traits::RangeType;
234 using RT_P =
typename LFSV_P::Traits::FiniteElementType::
235 Traits::LocalBasisType::Traits::RangeType;
238 const auto& lfsv_v_pfs =
child(lfsv,
_0);
239 const unsigned int vsize = lfsv_v_pfs.child(0).size();
240 const auto& lfsv_p =
child(lfsv,
_1);
241 const unsigned int psize = lfsv_p.size();
244 const int dim = EG::Geometry::mydimension;
247 const auto& cell = eg.entity();
250 auto geo = eg.geometry();
253 const int v_order = lfsv_v_pfs.child(0).finiteElement().localBasis().order();
254 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
255 const int qorder = 2*v_order + det_jac_order + superintegration_order;
258 std::vector<RT_V> phi(vsize);
259 std::vector<RT_P> psi(psize);
262 for (
const auto& ip : quadratureRule(geo,qorder))
264 lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
266 lfsv_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
269 const auto f1 = _p.f(cell,ip.position());
272 const auto factor = ip.weight() * geo.integrationElement(ip.position());
274 for(
int d=0; d<dim; ++d){
276 const auto& lfsv_v = lfsv_v_pfs.child(d);
278 for (
size_t i=0; i<vsize; i++)
281 r.accumulate(lfsv_v,i, -f1[d]*phi[i] * factor);
286 const auto g2 = _p.g2(eg,ip.position());
289 for (
size_t i=0; i<psize; i++)
291 r.accumulate(lfsv_p,i, g2 * psi[i] * factor);
299 template<
typename IG,
typename LFSV,
typename R>
300 void lambda_boundary (
const IG& ig,
const LFSV& lfsv, R& r)
const
304 using LFSV_V_PFS = TypeTree::Child<LFSV,_0>;
305 using LFSV_V = TypeTree::Child<LFSV_V_PFS,_0>;
306 using RT_V =
typename LFSV_V::Traits::FiniteElementType::
307 Traits::LocalBasisType::Traits::RangeType;
310 const auto& lfsv_v_pfs =
child(lfsv,
_0);
311 const unsigned int vsize = lfsv_v_pfs.child(0).size();
314 static const unsigned int dim = IG::Entity::dimension;
317 auto geo = ig.geometry();
320 auto geo_in_inside = ig.geometryInInside();
323 const int v_order = lfsv_v_pfs.child(0).finiteElement().localBasis().order();
324 const int det_jac_order = geo_in_inside.type().isSimplex() ? 0 : (dim-2);
325 const int jac_order = geo_in_inside.type().isSimplex() ? 0 : 1;
326 const int qorder = 2*v_order + det_jac_order + jac_order + superintegration_order;
329 std::vector<RT_V> phi(vsize);
332 for (
const auto& ip : quadratureRule(geo,qorder))
335 auto bctype = _p.bctype(ig,ip.position());
338 if (bctype != BC::StressNeumann)
342 auto local = geo_in_inside.global(ip.position());
345 lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(local,phi);
347 const auto factor = ip.weight() * geo.integrationElement(ip.position());
348 const auto normal = ig.unitOuterNormal(ip.position());
351 const auto neumann_stress = _p.j(ig,ip.position(),normal);
353 for(
unsigned int d=0; d<dim; ++d)
356 const auto& 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>
369 void jacobian_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv,
374 using LFSU_V_PFS = TypeTree::Child<LFSU,_0>;
375 using LFSU_V = TypeTree::Child<LFSU_V_PFS,_0>;
376 using LFSU_P = TypeTree::Child<LFSU,_1>;
377 using RF =
typename LFSU_V::Traits::FiniteElementType::
378 Traits::LocalBasisType::Traits::RangeFieldType;
379 using RT_V =
typename LFSU_V::Traits::FiniteElementType::
380 Traits::LocalBasisType::Traits::RangeType;
381 using JacobianType_V =
typename LFSU_V::Traits::FiniteElementType::
382 Traits::LocalBasisType::Traits::JacobianType;
383 using RT_P =
typename LFSU_P::Traits::FiniteElementType::
384 Traits::LocalBasisType::Traits::RangeType;
387 const auto& lfsu_v_pfs =
child(lfsu,
_0);
388 const unsigned int vsize = lfsu_v_pfs.child(0).size();
389 const auto& lfsu_p =
child(lfsu,
_1);
390 const unsigned int psize = lfsu_p.size();
393 const int dim = EG::Geometry::mydimension;
396 auto geo = eg.geometry();
399 const int v_order = lfsu_v_pfs.child(0).finiteElement().localBasis().order();
400 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
401 const int jac_order = geo.type().isSimplex() ? 0 : 1;
402 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
405 typename EG::Geometry::JacobianInverseTransposed jac;
406 std::vector<JacobianType_V> js(vsize);
407 std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
408 std::vector<RT_P> psi(psize);
409 std::vector<RT_V> phi(vsize);
414 for (
const auto& ip : quadratureRule(geo,qorder))
417 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(ip.position(),js);
420 jac = geo.jacobianInverseTransposed(ip.position());
421 for (
size_t i=0; i<vsize; i++)
424 jac.umv(js[i][0],gradphi[i]);
428 lfsu_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
432 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
434 for(
int d = 0; d < dim; ++d){
436 const auto& lfsv_v = lfsu_v_pfs.child(d);
437 for(
size_t l = 0; l < vsize; ++l)
438 vu[d] += x(lfsv_v,l) * phi[l];
443 const auto mu = _p.mu(eg,ip.position());
444 const auto rho = _p.rho(eg,ip.position());
446 const auto factor = ip.weight() * geo.integrationElement(ip.position());
448 for(
int d=0; d<dim; ++d){
450 const auto& lfsv_v = lfsu_v_pfs.child(d);
451 const auto& lfsu_v = lfsv_v;
456 for(
size_t l =0; l < vsize; ++l)
457 gradu_d.axpy(x(lfsv_v,l), gradphi[l]);
459 for (
size_t i=0; i<lfsv_v.size(); i++){
462 for (
size_t j=0; j<lfsv_v.size(); j++){
463 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[i] * gradphi[j]) * factor);
467 for(
int dd=0; dd<dim; ++dd){
468 const auto& lfsu_v = lfsu_v_pfs.child(dd);
469 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[j][d] * gradphi[i][dd]) * factor);
475 for (
size_t j=0; j<psize; j++)
476 mat.accumulate(lfsv_v,i,lfsu_p,j, - (gradphi[i][d] * psi[j]) * factor);
479 for(
int k =0; k < dim; ++k){
480 const auto& lfsu_v = lfsu_v_pfs.child(k);
482 const auto pre_factor = factor * rho * gradu_d[k] * phi[i];
484 for(
size_t j=0; j< lfsu_v.size(); ++j)
485 mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * phi[j]);
488 const auto 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);
505 const int superintegration_order;
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
sparsity pattern generator
Definition: pattern.hh:14
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
void setTime(R t_)
set time for subsequent evaluation
Definition: idefault.hh:104
Default flags for all local operators.
Definition: flags.hh:19
A local operator for the Navier-Stokes equations.
Definition: taylorhoodnavierstokes.hh:56
A few common exception classes.
Implements a vector constructed from a given type representing a field and a compile-time given size.
constexpr index_constant< 0 > _0
Compile time index with value 0.
Definition: indices.hh:52
constexpr index_constant< 1 > _1
Compile time index with value 1.
Definition: indices.hh:55
Namespace with predefined compile time indices for the range [0,19].
Definition: indices.hh:50
Dune namespace.
Definition: alignedallocator.hh:13
Definition: stokesparameter.hh:32
A unique label for each type of element that can occur in a grid.