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 TypeTree::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 TypeTree::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 auto geo = eg.geometry();
246 const int v_order = lfsv_v_pfs.child(0).finiteElement().localBasis().order();
247 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
248 const int qorder = 2*v_order + det_jac_order + superintegration_order;
251 std::vector<RT_V> phi(vsize);
252 std::vector<RT_P> psi(psize);
257 lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
259 lfsv_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
262 const auto f1 = _p.f(eg,ip.position());
265 const auto factor = ip.weight() * geo.integrationElement(ip.position());
267 for(
int d=0; d<
dim; ++d){
269 const auto& lfsv_v = lfsv_v_pfs.child(d);
271 for (
size_t i=0; i<vsize; i++)
274 r.accumulate(lfsv_v,i, -f1[d]*phi[i] * factor);
279 const auto g2 = _p.g2(eg,ip.position());
282 for (
size_t i=0; i<psize; i++)
284 r.accumulate(lfsv_p,i, g2 * psi[i] * factor);
292 template<
typename IG,
typename LFSV,
typename R>
296 using namespace TypeTree::Indices;
297 using LFSV_V_PFS = TypeTree::Child<LFSV,_0>;
298 using LFSV_V = TypeTree::Child<LFSV_V_PFS,_0>;
299 using RT_V =
typename LFSV_V::Traits::FiniteElementType::
300 Traits::LocalBasisType::Traits::RangeType;
303 const auto& lfsv_v_pfs = child(lfsv,_0);
304 const unsigned int vsize = lfsv_v_pfs.child(0).size();
307 static const unsigned int dim = IG::dimension;
310 auto geo = ig.geometry();
313 auto geo_in_inside = ig.geometryInInside();
316 const int v_order = lfsv_v_pfs.child(0).finiteElement().localBasis().order();
317 const int det_jac_order = geo_in_inside.type().isSimplex() ? 0 : (dim-2);
318 const int jac_order = geo_in_inside.type().isSimplex() ? 0 : 1;
319 const int qorder = 2*v_order + det_jac_order + jac_order + superintegration_order;
322 std::vector<RT_V> phi(vsize);
328 auto bctype = _p.bctype(ig,ip.position());
335 auto local = geo_in_inside.global(ip.position());
338 lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(local,phi);
340 const auto factor = ip.weight() * geo.integrationElement(ip.position());
341 const auto normal = ig.unitOuterNormal(ip.position());
344 const auto neumann_stress = _p.j(ig,ip.position(),normal);
346 for(
unsigned int d=0; d<
dim; ++d)
349 const auto& lfsv_v = lfsv_v_pfs.child(d);
351 for (
size_t i=0; i<vsize; i++)
353 r.accumulate(lfsv_v,i, neumann_stress[d] * phi[i] * factor);
361 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
366 using namespace TypeTree::Indices;
367 using LFSU_V_PFS = TypeTree::Child<LFSU,_0>;
368 using LFSU_V = TypeTree::Child<LFSU_V_PFS,_0>;
369 using LFSU_P = TypeTree::Child<LFSU,_1>;
370 using RF =
typename LFSU_V::Traits::FiniteElementType::
371 Traits::LocalBasisType::Traits::RangeFieldType;
372 using RT_V =
typename LFSU_V::Traits::FiniteElementType::
373 Traits::LocalBasisType::Traits::RangeType;
374 using JacobianType_V =
typename LFSU_V::Traits::FiniteElementType::
375 Traits::LocalBasisType::Traits::JacobianType;
376 using RT_P =
typename LFSU_P::Traits::FiniteElementType::
377 Traits::LocalBasisType::Traits::RangeType;
380 const auto& lfsu_v_pfs = child(lfsu,_0);
381 const unsigned int vsize = lfsu_v_pfs.child(0).size();
382 const auto& lfsu_p = child(lfsu,_1);
383 const unsigned int psize = lfsu_p.size();
386 const int dim = EG::Geometry::mydimension;
389 auto geo = eg.geometry();
392 const int v_order = lfsu_v_pfs.child(0).finiteElement().localBasis().order();
393 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
394 const int jac_order = geo.type().isSimplex() ? 0 : 1;
395 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
398 typename EG::Geometry::JacobianInverseTransposed jac;
399 std::vector<JacobianType_V> js(vsize);
400 std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
401 std::vector<RT_P> psi(psize);
402 std::vector<RT_V> phi(vsize);
403 Dune::FieldVector<RF,dim> vu(0.0);
404 Dune::FieldVector<RF,dim> gradu_d(0.0);
410 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(ip.position(),js);
413 jac = geo.jacobianInverseTransposed(ip.position());
414 for (
size_t i=0; i<vsize; i++)
417 jac.umv(js[i][0],gradphi[i]);
421 lfsu_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
425 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
427 for(
int d = 0; d <
dim; ++d){
429 const auto& lfsv_v = lfsu_v_pfs.child(d);
430 for(
size_t l = 0; l < vsize; ++l)
431 vu[d] += x(lfsv_v,l) * phi[l];
436 const auto mu = _p.mu(eg,ip.position());
437 const auto rho = _p.rho(eg,ip.position());
439 const auto factor = ip.weight() * geo.integrationElement(ip.position());
441 for(
int d=0; d<
dim; ++d){
443 const auto& lfsv_v = lfsu_v_pfs.child(d);
444 const auto& lfsu_v = lfsv_v;
449 for(
size_t l =0; l < vsize; ++l)
450 gradu_d.axpy(x(lfsv_v,l), gradphi[l]);
452 for (
size_t i=0; i<lfsv_v.size(); i++){
455 for (
size_t j=0; j<lfsv_v.size(); j++){
456 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[i] * gradphi[j]) * factor);
460 for(
int dd=0; dd<
dim; ++dd){
461 const auto& lfsu_v = lfsu_v_pfs.child(dd);
462 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[j][d] * gradphi[i][dd]) * factor);
468 for (
size_t j=0; j<psize; j++)
469 mat.accumulate(lfsv_v,i,lfsu_p,j, - (gradphi[i][d] * psi[j]) * factor);
472 for(
int k =0; k <
dim; ++k){
473 const auto& lfsu_v = lfsu_v_pfs.child(k);
475 const auto pre_factor = factor * rho * gradu_d[k] * phi[i];
477 for(
size_t j=0; j< lfsu_v.size(); ++j)
478 mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * phi[j]);
483 const auto pre_factor = factor * rho * phi[i];
484 for(
size_t j=0; j< lfsu_v.size(); ++j)
485 mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * (vu * gradphi[j]));
490 for (
size_t i=0; i<psize; i++){
491 for (
size_t j=0; j<lfsu_v.size(); j++)
492 mat.accumulate(lfsu_p,i,lfsu_v,j, - (gradphi[j][d] * psi[i]) * factor);
500 const int superintegration_order;
Definition: taylorhoodnavierstokes.hh:71
sparsity pattern generator
Definition: pattern.hh:13
static const bool navier
Definition: taylorhoodnavierstokes.hh:67
const IG & ig
Definition: constraints.hh:148
static const int dim
Definition: adaptivity.hh:83
Definition: taylorhoodnavierstokes.hh:76
Definition: adaptivity.hh:27
Definition: stokesparameter.hh:32
TaylorHoodNavierStokes(const PhysicalParameters &p, int superintegration_order_=0)
Definition: taylorhoodnavierstokes.hh:80
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:88
Default flags for all local operators.
Definition: flags.hh:18
void jacobian_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, M &mat) const
Definition: taylorhoodnavierstokes.hh:362
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
static const bool full_tensor
Definition: taylorhoodnavierstokes.hh:68
A local operator for the Navier-Stokes equations.
Definition: taylorhoodnavierstokes.hh:58
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
Definition: stokesparameter.hh:36
P PhysicalParameters
Definition: taylorhoodnavierstokes.hh:78
Definition: taylorhoodnavierstokes.hh:74
const P & p
Definition: constraints.hh:147
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:293
Definition: taylorhoodnavierstokes.hh:75
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:221