350 libmesh_assert_equal_to (system_name,
"Navier-Stokes");
364 const unsigned int u_var = navier_stokes_system.variable_number (
"vel_x");
365 const unsigned int v_var = navier_stokes_system.variable_number (
"vel_y");
366 const unsigned int p_var = navier_stokes_system.variable_number (
"p");
370 FEType fe_vel_type = navier_stokes_system.variable_type(u_var);
373 FEType fe_pres_type = navier_stokes_system.variable_type(p_var);
377 std::unique_ptr<FEBase> fe_vel (FEBase::build(
dim, fe_vel_type));
381 std::unique_ptr<FEBase> fe_pres (FEBase::build(
dim, fe_pres_type));
388 fe_vel->attach_quadrature_rule (&qrule);
389 fe_pres->attach_quadrature_rule (&qrule);
395 const std::vector<Real> & JxW = fe_vel->get_JxW();
398 const std::vector<std::vector<Real>> & phi = fe_vel->get_phi();
402 const std::vector<std::vector<RealGradient>> & dphi = fe_vel->get_dphi();
406 const std::vector<std::vector<Real>> & psi = fe_pres->get_phi();
415 const DofMap & dof_map = navier_stokes_system.get_dof_map();
425 Kuu(Ke), Kuv(Ke), Kup(Ke),
426 Kvu(Ke), Kvv(Ke), Kvp(Ke),
427 Kpu(Ke), Kpv(Ke), Kpp(Ke);
435 std::reference_wrapper<DenseSubVector<Number>> F[2] = {Fu, Fv};
438 std::reference_wrapper<DenseSubMatrix<Number>> K[2][2] = {{Kuu, Kuv}, {Kvu, Kvv}};
441 std::reference_wrapper<DenseSubMatrix<Number>>
B[2] = {Kup, Kvp};
442 std::reference_wrapper<DenseSubMatrix<Number>> BT[2] = {Kpu, Kpv};
447 std::vector<dof_id_type> dof_indices;
448 std::vector<dof_id_type> dof_indices_u;
449 std::vector<dof_id_type> dof_indices_v;
450 std::vector<dof_id_type> dof_indices_p;
463 const Real theta = 1.;
471 const bool pin_pressure = es.
parameters.
get<
bool>(
"pin_pressure");
481 for (
const auto & elem :
mesh.active_local_element_ptr_range())
492 const unsigned int n_dofs = dof_indices.size();
493 const unsigned int n_u_dofs = dof_indices_u.size();
494 const unsigned int n_v_dofs = dof_indices_v.size();
495 const unsigned int n_p_dofs = dof_indices_p.size();
501 fe_vel->reinit (elem);
502 fe_pres->reinit (elem);
510 Ke.
resize (n_dofs, n_dofs);
526 Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
527 Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
528 Kup.reposition (u_var*n_u_dofs, p_var*n_u_dofs, n_u_dofs, n_p_dofs);
530 Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
531 Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
532 Kvp.reposition (v_var*n_v_dofs, p_var*n_v_dofs, n_v_dofs, n_p_dofs);
534 Kpu.reposition (p_var*n_u_dofs, u_var*n_u_dofs, n_p_dofs, n_u_dofs);
535 Kpv.reposition (p_var*n_u_dofs, v_var*n_u_dofs, n_p_dofs, n_v_dofs);
536 Kpp.reposition (p_var*n_u_dofs, p_var*n_u_dofs, n_p_dofs, n_p_dofs);
538 Fu.reposition (u_var*n_u_dofs, n_u_dofs);
539 Fv.reposition (v_var*n_u_dofs, n_v_dofs);
540 Fp.reposition (p_var*n_u_dofs, n_p_dofs);
548 for (
unsigned int qp=0; qp<qrule.n_points(); qp++)
556 std::array<Gradient, 2> grad_uv{};
559 std::array<Gradient, 2> grad_uv_old{};
563 for (
unsigned int l=0; l<n_u_dofs; l++)
566 U_old(0) += phi[l][qp]*navier_stokes_system.
old_solution (dof_indices_u[l]);
567 U_old(1) += phi[l][qp]*navier_stokes_system.
old_solution (dof_indices_v[l]);
568 grad_uv_old[0].add_scaled (dphi[l][qp],navier_stokes_system.
old_solution (dof_indices_u[l]));
569 grad_uv_old[1].add_scaled (dphi[l][qp],navier_stokes_system.
old_solution (dof_indices_v[l]));
572 U(0) += phi[l][qp]*navier_stokes_system.current_solution (dof_indices_u[l]);
573 U(1) += phi[l][qp]*navier_stokes_system.current_solution (dof_indices_v[l]);
574 grad_uv[0].add_scaled (dphi[l][qp],navier_stokes_system.current_solution (dof_indices_u[l]));
575 grad_uv[1].add_scaled (dphi[l][qp],navier_stokes_system.current_solution (dof_indices_v[l]));
579 for (
unsigned int l=0; l<n_p_dofs; l++)
580 p_old += psi[l][qp]*navier_stokes_system.
old_solution (dof_indices_p[l]);
585 for (
unsigned int i=0; i<n_u_dofs; i++)
587 for (
unsigned int k=0; k<2; ++k)
589 (U_old(k) * phi[i][qp] -
590 (1.-theta) * dt * (U_old * grad_uv_old[k]) * phi[i][qp] +
591 (1.-theta) * dt * p_old * dphi[i][qp](k) -
592 (1.-theta) * dt * nu * (grad_uv_old[k] * dphi[i][qp]) +
593 theta * dt * (U * grad_uv[k]) * phi[i][qp]);
596 for (
unsigned int j=0; j<n_u_dofs; j++)
597 for (
unsigned int k=0; k<2; ++k)
598 for (
unsigned int l=0; l<2; ++l)
602 K[k][k](i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] +
603 theta*dt*nu*(dphi[i][qp]*dphi[j][qp]) +
604 theta*dt*(U*dphi[j][qp])*phi[i][qp]);
607 K[k][l](i,j) += JxW[qp] * theta * dt * grad_uv[k](l) * phi[i][qp] * phi[j][qp];
611 for (
unsigned int j=0; j<n_p_dofs; j++)
612 for (
unsigned int k=0; k<2; ++k)
613 B[k](i,j) += JxW[qp] * -theta * dt * psi[j][qp] * dphi[i][qp](k);
620 for (
unsigned int i=0; i<n_p_dofs; i++)
621 for (
unsigned int j=0; j<n_u_dofs; j++)
622 for (
unsigned int k=0; k<2; ++k)
623 BT[k](i,j) += JxW[qp] * psi[i][qp] * dphi[j][qp](k);
637 const Real penalty = 1.e10;
638 const unsigned int pressure_node = 0;
639 const Real p_value = 0.0;
640 for (
auto c : elem->node_index_range())
641 if (elem->node_id(c) == pressure_node)
644 Fp(c) += penalty*p_value;
657 matrix.add_matrix (Ke, dof_indices);
658 navier_stokes_system.rhs->add_vector (Fe, dof_indices);
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
void resize(const unsigned int n)
Resize the vector.
Manages storage and variables for transient systems.
Order default_quadrature_order() const
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
const T_sys & get_system(std::string_view name) const
This is the MeshBase class.
Defines a dense subvector for use in finite element computations.
This class handles the numbering of degrees of freedom on a mesh.
void libmesh_ignore(const Args &...)
const T & get(std::string_view) const
Defines a dense submatrix for use in Finite Element-type computations.
Number old_solution(const dof_id_type global_dof_number) const
void heterogenously_constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true, int qoi_index=-1) const
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const MeshBase & get_mesh() const
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
This class implements specific orders of Gauss quadrature.
unsigned int mesh_dimension() const
Parameters parameters
Data structure holding arbitrary parameters.