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> 70 static const bool navier = P::assemble_navier;
86 , superintegration_order(superintegration_order_)
97 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename R>
98 void alpha_volume (
const EG& eg,
const LFSU& lfsu,
const X& x,
const LFSV& lfsv, R& r)
const 101 using namespace Indices;
102 using LFSU_V_PFS = TypeTree::Child<LFSU,_0>;
103 using LFSU_V = TypeTree::Child<LFSU_V_PFS,_0>;
104 using LFSU_P = TypeTree::Child<LFSU,_1>;
105 using RF =
typename LFSU_V::Traits::FiniteElementType::
106 Traits::LocalBasisType::Traits::RangeFieldType;
107 using RT_V =
typename LFSU_V::Traits::FiniteElementType::
108 Traits::LocalBasisType::Traits::RangeType;
109 using JacobianType_V =
typename LFSU_V::Traits::FiniteElementType::
110 Traits::LocalBasisType::Traits::JacobianType;
111 using RT_P =
typename LFSU_P::Traits::FiniteElementType::
112 Traits::LocalBasisType::Traits::RangeType;
115 const auto& lfsu_v_pfs = child(lfsu,_0);
116 const unsigned int vsize = lfsu_v_pfs.child(0).size();
117 const auto& lfsu_p = child(lfsu,_1);
118 const unsigned int psize = lfsu_p.size();
121 const int dim = EG::Geometry::mydimension;
124 auto geo = eg.geometry();
127 const int v_order = lfsu_v_pfs.child(0).finiteElement().localBasis().order();
128 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
129 const int jac_order = geo.type().isSimplex() ? 0 : 1;
130 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
133 typename EG::Geometry::JacobianInverseTransposed jac;
134 std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
135 std::vector<RT_P> psi(psize);
136 Dune::FieldVector<RF,dim> vu(0.0);
137 std::vector<RT_V> phi(vsize);
138 Dune::FieldMatrix<RF,dim,dim> jacu(0.0);
144 std::vector<JacobianType_V> js(vsize);
145 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(ip.position(),js);
148 jac = geo.jacobianInverseTransposed(ip.position());
149 for (
size_t i=0; i<vsize; i++)
152 jac.umv(js[i][0],gradphi[i]);
156 lfsu_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
161 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
163 for(
int d=0; d<
dim; ++d)
166 const auto& lfsu_v = lfsu_v_pfs.child(d);
167 for (
size_t i=0; i<lfsu_v.size(); i++)
168 vu[d] += x(lfsu_v,i) * phi[i];
173 for(
int d=0; d<
dim; ++d){
175 const auto& lfsu_v = lfsu_v_pfs.child(d);
176 for (
size_t i=0; i<lfsu_v.size(); i++)
177 jacu[d].axpy(x(lfsu_v,i),gradphi[i]);
182 for (
size_t i=0; i<lfsu_p.size(); i++)
183 func_p += x(lfsu_p,i) * psi[i];
186 const auto mu = _p.mu(eg,ip.position());
187 const auto rho = _p.rho(eg,ip.position());
190 const auto factor = ip.weight() * geo.integrationElement(ip.position());
192 for(
int d=0; d<
dim; ++d){
194 const auto& lfsu_v = lfsu_v_pfs.child(d);
197 const auto u_nabla_u = vu * jacu[d];
199 for (
size_t i=0; i<vsize; i++){
202 r.accumulate(lfsu_v,i, mu * (jacu[d] * gradphi[i]) * factor);
206 for(
int dd=0; dd<
dim; ++dd)
207 r.accumulate(lfsu_v,i, mu * (jacu[dd][d] * gradphi[i][dd]) * factor);
210 r.accumulate(lfsu_v,i,- (func_p * gradphi[i][d]) * factor);
214 r.accumulate(lfsu_v,i, rho * u_nabla_u * phi[i] * factor);
220 for (
size_t i=0; i<psize; i++)
221 for(
int d=0; d<
dim; ++d)
223 r.accumulate(lfsu_p,i, -1.0 * jacu[d][d] * psi[i] * factor);
230 template<
typename EG,
typename LFSV,
typename R>
234 using namespace Indices;
235 using LFSV_V_PFS = TypeTree::Child<LFSV,_0>;
236 using LFSV_V = TypeTree::Child<LFSV_V_PFS,_0>;
237 using LFSV_P = TypeTree::Child<LFSV,_1>;
238 using RT_V =
typename LFSV_V::Traits::FiniteElementType::
239 Traits::LocalBasisType::Traits::RangeType;
240 using RT_P =
typename LFSV_P::Traits::FiniteElementType::
241 Traits::LocalBasisType::Traits::RangeType;
244 const auto& lfsv_v_pfs = child(lfsv,_0);
245 const unsigned int vsize = lfsv_v_pfs.child(0).size();
246 const auto& lfsv_p = child(lfsv,_1);
247 const unsigned int psize = lfsv_p.size();
250 const int dim = EG::Geometry::mydimension;
253 const auto& cell = eg.entity();
256 auto geo = eg.geometry();
259 const int v_order = lfsv_v_pfs.child(0).finiteElement().localBasis().order();
260 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
261 const int qorder = 2*v_order + det_jac_order + superintegration_order;
264 std::vector<RT_V> phi(vsize);
265 std::vector<RT_P> psi(psize);
270 lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
272 lfsv_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
275 const auto f1 = _p.f(cell,ip.position());
278 const auto factor = ip.weight() * geo.integrationElement(ip.position());
280 for(
int d=0; d<
dim; ++d){
282 const auto& lfsv_v = lfsv_v_pfs.child(d);
284 for (
size_t i=0; i<vsize; i++)
287 r.accumulate(lfsv_v,i, -f1[d]*phi[i] * factor);
292 const auto g2 = _p.g2(eg,ip.position());
295 for (
size_t i=0; i<psize; i++)
297 r.accumulate(lfsv_p,i, g2 * psi[i] * factor);
305 template<
typename IG,
typename LFSV,
typename R>
309 using namespace Indices;
310 using LFSV_V_PFS = TypeTree::Child<LFSV,_0>;
311 using LFSV_V = TypeTree::Child<LFSV_V_PFS,_0>;
312 using RT_V =
typename LFSV_V::Traits::FiniteElementType::
313 Traits::LocalBasisType::Traits::RangeType;
316 const auto& lfsv_v_pfs = child(lfsv,_0);
317 const unsigned int vsize = lfsv_v_pfs.child(0).size();
320 static const unsigned int dim = IG::dimension;
323 auto geo = ig.geometry();
326 auto geo_in_inside = ig.geometryInInside();
329 const int v_order = lfsv_v_pfs.child(0).finiteElement().localBasis().order();
330 const int det_jac_order = geo_in_inside.type().isSimplex() ? 0 : (dim-2);
331 const int jac_order = geo_in_inside.type().isSimplex() ? 0 : 1;
332 const int qorder = 2*v_order + det_jac_order + jac_order + superintegration_order;
335 std::vector<RT_V> phi(vsize);
341 auto bctype = _p.bctype(ig,ip.position());
348 auto local = geo_in_inside.global(ip.position());
351 lfsv_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(local,phi);
353 const auto factor = ip.weight() * geo.integrationElement(ip.position());
354 const auto normal = ig.unitOuterNormal(ip.position());
357 const auto neumann_stress = _p.j(ig,ip.position(),normal);
359 for(
unsigned int d=0; d<
dim; ++d)
362 const auto& lfsv_v = lfsv_v_pfs.child(d);
364 for (
size_t i=0; i<vsize; i++)
366 r.accumulate(lfsv_v,i, neumann_stress[d] * phi[i] * factor);
374 template<
typename EG,
typename LFSU,
typename X,
typename LFSV,
typename M>
379 using namespace Indices;
380 using LFSU_V_PFS = TypeTree::Child<LFSU,_0>;
381 using LFSU_V = TypeTree::Child<LFSU_V_PFS,_0>;
382 using LFSU_P = TypeTree::Child<LFSU,_1>;
383 using RF =
typename LFSU_V::Traits::FiniteElementType::
384 Traits::LocalBasisType::Traits::RangeFieldType;
385 using RT_V =
typename LFSU_V::Traits::FiniteElementType::
386 Traits::LocalBasisType::Traits::RangeType;
387 using JacobianType_V =
typename LFSU_V::Traits::FiniteElementType::
388 Traits::LocalBasisType::Traits::JacobianType;
389 using RT_P =
typename LFSU_P::Traits::FiniteElementType::
390 Traits::LocalBasisType::Traits::RangeType;
393 const auto& lfsu_v_pfs = child(lfsu,_0);
394 const unsigned int vsize = lfsu_v_pfs.child(0).size();
395 const auto& lfsu_p = child(lfsu,_1);
396 const unsigned int psize = lfsu_p.size();
399 const int dim = EG::Geometry::mydimension;
402 auto geo = eg.geometry();
405 const int v_order = lfsu_v_pfs.child(0).finiteElement().localBasis().order();
406 const int det_jac_order = geo.type().isSimplex() ? 0 : (dim-1);
407 const int jac_order = geo.type().isSimplex() ? 0 : 1;
408 const int qorder = 3*v_order - 1 + jac_order + det_jac_order + superintegration_order;
411 typename EG::Geometry::JacobianInverseTransposed jac;
412 std::vector<JacobianType_V> js(vsize);
413 std::vector<Dune::FieldVector<RF,dim> > gradphi(vsize);
414 std::vector<RT_P> psi(psize);
415 std::vector<RT_V> phi(vsize);
416 Dune::FieldVector<RF,dim> vu(0.0);
417 Dune::FieldVector<RF,dim> gradu_d(0.0);
423 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateJacobian(ip.position(),js);
426 jac = geo.jacobianInverseTransposed(ip.position());
427 for (
size_t i=0; i<vsize; i++)
430 jac.umv(js[i][0],gradphi[i]);
434 lfsu_p.finiteElement().localBasis().evaluateFunction(ip.position(),psi);
438 lfsu_v_pfs.child(0).finiteElement().localBasis().evaluateFunction(ip.position(),phi);
440 for(
int d = 0; d <
dim; ++d){
442 const auto& lfsv_v = lfsu_v_pfs.child(d);
443 for(
size_t l = 0; l < vsize; ++l)
444 vu[d] += x(lfsv_v,l) * phi[l];
449 const auto mu = _p.mu(eg,ip.position());
450 const auto rho = _p.rho(eg,ip.position());
452 const auto factor = ip.weight() * geo.integrationElement(ip.position());
454 for(
int d=0; d<
dim; ++d){
456 const auto& lfsv_v = lfsu_v_pfs.child(d);
457 const auto& lfsu_v = lfsv_v;
462 for(
size_t l =0; l < vsize; ++l)
463 gradu_d.axpy(x(lfsv_v,l), gradphi[l]);
465 for (
size_t i=0; i<lfsv_v.size(); i++){
468 for (
size_t j=0; j<lfsv_v.size(); j++){
469 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[i] * gradphi[j]) * factor);
473 for(
int dd=0; dd<
dim; ++dd){
474 const auto& lfsu_v = lfsu_v_pfs.child(dd);
475 mat.accumulate(lfsv_v,i,lfsu_v,j, mu * (gradphi[j][d] * gradphi[i][dd]) * factor);
481 for (
size_t j=0; j<psize; j++)
482 mat.accumulate(lfsv_v,i,lfsu_p,j, - (gradphi[i][d] * psi[j]) * factor);
485 for(
int k =0; k <
dim; ++k){
486 const auto& lfsu_v = lfsu_v_pfs.child(k);
488 const auto pre_factor = factor * rho * gradu_d[k] * phi[i];
490 for(
size_t j=0; j< lfsu_v.size(); ++j)
491 mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * phi[j]);
496 const auto pre_factor = factor * rho * phi[i];
497 for(
size_t j=0; j< lfsu_v.size(); ++j)
498 mat.accumulate(lfsv_v,i,lfsu_v,j, pre_factor * (vu * gradphi[j]));
503 for (
size_t i=0; i<psize; i++){
504 for (
size_t j=0; j<lfsu_v.size(); j++)
505 mat.accumulate(lfsu_p,i,lfsu_v,j, - (gradphi[j][d] * psi[i]) * factor);
513 const int superintegration_order;
520 #endif // DUNE_PDELAB_LOCALOPERATOR_TAYLORHOODNAVIERSTOKES_HH typename InstatBase::RealType Real
Definition: taylorhoodnavierstokes.hh:68
const IG & ig
Definition: constraints.hh:148
static const int dim
Definition: adaptivity.hh:83
Definition: taylorhoodnavierstokes.hh:78
void lambda_boundary(const IG &ig, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:306
Definition: stokesparameter.hh:36
Definition: taylorhoodnavierstokes.hh:77
static const bool full_tensor
Definition: taylorhoodnavierstokes.hh:71
void setTime(R t_)
set time for subsequent evaluation
Definition: idefault.hh:104
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:98
const P & p
Definition: constraints.hh:147
P PhysicalParameters
Definition: taylorhoodnavierstokes.hh:81
Definition: taylorhoodnavierstokes.hh:74
TaylorHoodNavierStokes(PhysicalParameters &p, int superintegration_order_=0)
Definition: taylorhoodnavierstokes.hh:83
void setTime(Real t)
Definition: taylorhoodnavierstokes.hh:90
Definition: stokesparameter.hh:32
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:375
A local operator for the Navier-Stokes equations.
Definition: taylorhoodnavierstokes.hh:58
void lambda_volume(const EG &eg, const LFSV &lfsv, R &r) const
Definition: taylorhoodnavierstokes.hh:231
Default flags for all local operators.
Definition: flags.hh:18
static const bool navier
Definition: taylorhoodnavierstokes.hh:70
Definition: taylorhoodnavierstokes.hh:79
R RealType
Definition: idefault.hh:92