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> 67 static const bool navier = P::assemble_navier;
83 , superintegration_order(superintegration_order_)
87 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
88 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 91 using namespace Indices;
92 using LFSU_V_PFS = TypeTree::Child<LFSU,_0>;
93 using LFSU_V = TypeTree::Child<LFSU_V_PFS,_0>;
94 using LFSU_P = TypeTree::Child<LFSU,_1>;
95 using RF =
typename LFSU_V::Traits::FiniteElementType::
96 Traits::LocalBasisType::Traits::RangeFieldType;
97 using RT_V =
typename LFSU_V::Traits::FiniteElementType::
98 Traits::LocalBasisType::Traits::RangeType;
99 using JacobianType_V =
typename LFSU_V::Traits::FiniteElementType::
100 Traits::LocalBasisType::Traits::JacobianType;
101 using RT_P =
typename LFSU_P::Traits::FiniteElementType::
102 Traits::LocalBasisType::Traits::RangeType;
105 const auto& lfsu_v_pfs = child(lfsu,_0);
106 const unsigned int vsize = lfsu_v_pfs.child(0).size();
107 const auto& lfsu_p = child(lfsu,_1);
108 const unsigned int psize = lfsu_p.size();
111 const int dim = EG::Geometry::mydimension;
114 auto geo = eg.geometry();
117 const int v_order = lfsu_v_pfs.child(0).finiteElement().localBasis().order();
118 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
119 const int jac_order = geo.type().isSimplex() ? 0 : 1;
120 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
123 typename EG::Geometry::JacobianInverseTransposed jac;
124 std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
125 std::vector<RT_P> psi(psize);
126 Dune::FieldVector<RF,dim> vu(0.0);
127 std::vector<RT_V> phi(vsize);
128 Dune::FieldMatrix<RF,dim,dim> jacu(0.0);
134 std::vector<JacobianType_V> js(vsize);
135 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(ip.position(),js);
138 jac = geo.jacobianInverseTransposed(ip.position());
139 for (
size_t i=0; i<vsize; i++)
142 jac.umv(js[i][0],gradphi[i]);
146 lfsu_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
151 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
153 for(
int d=0; d<
dim; ++d)
156 const auto& lfsu_v = lfsu_v_pfs.child(d);
157 for (
size_t i=0; i<lfsu_v.size(); i++)
158 vu[d] += x(lfsu_v,i) * phi[i];
163 for(
int d=0; d<
dim; ++d){
165 const auto& lfsu_v = lfsu_v_pfs.child(d);
166 for (
size_t i=0; i<lfsu_v.size(); i++)
167 jacu[d].axpy(x(lfsu_v,i),gradphi[i]);
172 for (
size_t i=0; i<lfsu_p.size(); i++)
173 func_p += x(lfsu_p,i) * psi[i];
176 const auto mu = _p.mu(eg,ip.position());
177 const auto rho = _p.rho(eg,ip.position());
180 const auto factor = ip.weight() * geo.integrationElement(ip.position());
182 for(
int d=0; d<
dim; ++d){
184 const auto& lfsu_v = lfsu_v_pfs.child(d);
187 const auto u_nabla_u = vu * jacu[d];
189 for (
size_t i=0; i<vsize; i++){
192 r.accumulate(lfsu_v,i, mu * (jacu[d] * gradphi[i]) * factor);
196 for(
int dd=0; dd<
dim; ++dd)
197 r.accumulate(lfsu_v,i, mu * (jacu[dd][d] * gradphi[i][dd]) * factor);
200 r.accumulate(lfsu_v,i,- (func_p * gradphi[i][d]) * factor);
204 r.accumulate(lfsu_v,i, rho * u_nabla_u * phi[i] * factor);
210 for (
size_t i=0; i<psize; i++)
211 for(
int d=0; d<
dim; ++d)
213 r.accumulate(lfsu_p,i, -1.0 * jacu[d][d] * psi[i] * factor);
220 template<
typename EG,
typename LFSV,
typename R>
224 using namespace Indices;
225 using LFSV_V_PFS = TypeTree::Child<LFSV,_0>;
226 using LFSV_V = TypeTree::Child<LFSV_V_PFS,_0>;
227 using LFSV_P = TypeTree::Child<LFSV,_1>;
228 using RT_V =
typename LFSV_V::Traits::FiniteElementType::
229 Traits::LocalBasisType::Traits::RangeType;
230 using RT_P =
typename LFSV_P::Traits::FiniteElementType::
231 Traits::LocalBasisType::Traits::RangeType;
234 const auto& lfsv_v_pfs = child(lfsv,_0);
235 const unsigned int vsize = lfsv_v_pfs.child(0).size();
236 const auto& lfsv_p = child(lfsv,_1);
237 const unsigned int psize = lfsv_p.size();
240 const int dim = EG::Geometry::mydimension;
243 const auto& cell = eg.entity();
246 auto geo = eg.geometry();
249 const int v_order = lfsv_v_pfs.child(0).finiteElement().localBasis().order();
250 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
251 const int qorder = 2*v_order + det_jac_order + superintegration_order;
254 std::vector<RT_V> phi(vsize);
255 std::vector<RT_P> psi(psize);
260 lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
262 lfsv_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
265 const auto f1 = _p.f(cell,ip.position());
268 const auto factor = ip.weight() * geo.integrationElement(ip.position());
270 for(
int d=0; d<
dim; ++d){
272 const auto& lfsv_v = lfsv_v_pfs.child(d);
274 for (
size_t i=0; i<vsize; i++)
277 r.accumulate(lfsv_v,i, -f1[d]*phi[i] * factor);
282 const auto g2 = _p.g2(eg,ip.position());
285 for (
size_t i=0; i<psize; i++)
287 r.accumulate(lfsv_p,i, g2 * psi[i] * factor);
295 template<
typename IG,
typename LFSV,
typename R>
299 using namespace Indices;
300 using LFSV_V_PFS = TypeTree::Child<LFSV,_0>;
301 using LFSV_V = TypeTree::Child<LFSV_V_PFS,_0>;
302 using RT_V =
typename LFSV_V::Traits::FiniteElementType::
303 Traits::LocalBasisType::Traits::RangeType;
306 const auto& lfsv_v_pfs = child(lfsv,_0);
307 const unsigned int vsize = lfsv_v_pfs.child(0).size();
310 static const unsigned int dim = IG::dimension;
313 auto geo = ig.geometry();
316 auto geo_in_inside = ig.geometryInInside();
319 const int v_order = lfsv_v_pfs.child(0).finiteElement().localBasis().order();
320 const int det_jac_order = geo_in_inside.type().isSimplex() ? 0 : (dim-2);
321 const int jac_order = geo_in_inside.type().isSimplex() ? 0 : 1;
322 const int qorder = 2*v_order + det_jac_order + jac_order + superintegration_order;
325 std::vector<RT_V> phi(vsize);
331 auto bctype = _p.bctype(ig,ip.position());
338 auto local = geo_in_inside.global(ip.position());
341 lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(local,phi);
343 const auto factor = ip.weight() * geo.integrationElement(ip.position());
344 const auto normal = ig.unitOuterNormal(ip.position());
347 const auto neumann_stress = _p.j(ig,ip.position(),normal);
349 for(
unsigned int d=0; d<
dim; ++d)
352 const auto& 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 using namespace Indices;
370 using LFSU_V_PFS = TypeTree::Child<LFSU,_0>;
371 using LFSU_V = TypeTree::Child<LFSU_V_PFS,_0>;
372 using LFSU_P = TypeTree::Child<LFSU,_1>;
373 using RF =
typename LFSU_V::Traits::FiniteElementType::
374 Traits::LocalBasisType::Traits::RangeFieldType;
375 using RT_V =
typename LFSU_V::Traits::FiniteElementType::
376 Traits::LocalBasisType::Traits::RangeType;
377 using JacobianType_V =
typename LFSU_V::Traits::FiniteElementType::
378 Traits::LocalBasisType::Traits::JacobianType;
379 using RT_P =
typename LFSU_P::Traits::FiniteElementType::
380 Traits::LocalBasisType::Traits::RangeType;
383 const auto& lfsu_v_pfs = child(lfsu,_0);
384 const unsigned int vsize = lfsu_v_pfs.child(0).size();
385 const auto& lfsu_p = child(lfsu,_1);
386 const unsigned int psize = lfsu_p.size();
389 const int dim = EG::Geometry::mydimension;
392 auto geo = eg.geometry();
395 const int v_order = lfsu_v_pfs.child(0).finiteElement().localBasis().order();
396 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
397 const int jac_order = geo.type().isSimplex() ? 0 : 1;
398 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
401 typename EG::Geometry::JacobianInverseTransposed jac;
402 std::vector<JacobianType_V> js(vsize);
403 std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
404 std::vector<RT_P> psi(psize);
405 std::vector<RT_V> phi(vsize);
406 Dune::FieldVector<RF,dim> vu(0.0);
407 Dune::FieldVector<RF,dim> gradu_d(0.0);
413 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(ip.position(),js);
416 jac = geo.jacobianInverseTransposed(ip.position());
417 for (
size_t i=0; i<vsize; i++)
420 jac.umv(js[i][0],gradphi[i]);
424 lfsu_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
428 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
430 for(
int d = 0; d <
dim; ++d){
432 const auto& lfsv_v = lfsu_v_pfs.child(d);
433 for(
size_t l = 0; l < vsize; ++l)
434 vu[d] += x(lfsv_v,l) * phi[l];
439 const auto mu = _p.mu(eg,ip.position());
440 const auto rho = _p.rho(eg,ip.position());
442 const auto factor = ip.weight() * geo.integrationElement(ip.position());
444 for(
int d=0; d<
dim; ++d){
446 const auto& lfsv_v = lfsu_v_pfs.child(d);
447 const auto& lfsu_v = lfsv_v;
452 for(
size_t l =0; l < vsize; ++l)
453 gradu_d.axpy(x(lfsv_v,l), gradphi[l]);
455 for (
size_t i=0; i<lfsv_v.size(); i++){
458 for (
size_t j=0; j<lfsv_v.size(); j++){
459 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[i] * gradphi[j]) * factor);
463 for(
int dd=0; dd<
dim; ++dd){
464 const auto& lfsu_v = lfsu_v_pfs.child(dd);
465 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[j][d] * gradphi[i][dd]) * factor);
471 for (
size_t j=0; j<psize; j++)
472 mat.accumulate(lfsv_v,i,lfsu_p,j, - (gradphi[i][d] * psi[j]) * factor);
475 for(
int k =0; k <
dim; ++k){
476 const auto& lfsu_v = lfsu_v_pfs.child(k);
478 const auto pre_factor = factor * rho * gradu_d[k] * phi[i];
480 for(
size_t j=0; j< lfsu_v.size(); ++j)
481 mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * phi[j]);
486 const auto pre_factor = factor * rho * phi[i];
487 for(
size_t j=0; j< lfsu_v.size(); ++j)
488 mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * (vu * gradphi[j]));
493 for (
size_t i=0; i<psize; i++){
494 for (
size_t j=0; j<lfsu_v.size(); j++)
495 mat.accumulate(lfsu_p,i,lfsu_v,j, - (gradphi[j][d] * psi[i]) * factor);
503 const int superintegration_order;
510 #endif // DUNE_PDELAB_LOCALOPERATOR_TAYLORHOODNAVIERSTOKES_HH static const int dim
Definition: adaptivity.hh:83
const IG & ig
Definition: constraints.hh:148
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:88
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Definition: stokesparameter.hh:32
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
Definition: taylorhoodnavierstokes.hh:71
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:221
sparsity pattern generator
Definition: pattern.hh:13
QuadratureRuleWrapper< QuadratureRule< typename Geometry::ctype, Geometry::mydimension > > quadratureRule(const Geometry &geo, std::size_t order, QuadratureType::Enum quadrature_type=QuadratureType::GaussLegendre)
Returns a quadrature rule for the given geometry.
Definition: quadraturerules.hh:117
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: taylorhoodnavierstokes.hh:365
A local operator for the Navier-Stokes equations.
Definition: taylorhoodnavierstokes.hh:58
Definition: taylorhoodnavierstokes.hh:75
Default flags for all local operators.
Definition: flags.hh:18
Definition: stokesparameter.hh:36
Definition: taylorhoodnavierstokes.hh:76
static const bool full_tensor
Definition: taylorhoodnavierstokes.hh:68
P PhysicalParameters
Definition: taylorhoodnavierstokes.hh:78
TaylorHoodNavierStokes(const PhysicalParameters &p, int superintegration_order_=0)
Definition: taylorhoodnavierstokes.hh:80
const P & p
Definition: constraints.hh:147
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:296
Definition: taylorhoodnavierstokes.hh:74
static const bool navier
Definition: taylorhoodnavierstokes.hh:67