344 libmesh_assert_equal_to (system_name,
"Navier-Stokes");
358 const unsigned int u_var = navier_stokes_system.variable_number (
"vel_x");
359 const unsigned int v_var = navier_stokes_system.variable_number (
"vel_y");
360 const unsigned int p_var = navier_stokes_system.variable_number (
"p");
361 const unsigned int alpha_var = navier_stokes_system.variable_number (
"alpha");
365 FEType fe_vel_type = navier_stokes_system.variable_type(u_var);
368 FEType fe_pres_type = navier_stokes_system.variable_type(p_var);
372 std::unique_ptr<FEBase> fe_vel (FEBase::build(
dim, fe_vel_type));
376 std::unique_ptr<FEBase> fe_pres (FEBase::build(
dim, fe_pres_type));
383 fe_vel->attach_quadrature_rule (&qrule);
384 fe_pres->attach_quadrature_rule (&qrule);
390 const std::vector<Real> & JxW = fe_vel->get_JxW();
393 const std::vector<std::vector<Real>> & phi = fe_vel->get_phi();
397 const std::vector<std::vector<RealGradient>> & dphi = fe_vel->get_dphi();
401 const std::vector<std::vector<Real>> & psi = fe_pres->get_phi();
410 const DofMap & dof_map = navier_stokes_system.get_dof_map();
420 Kuu(Ke), Kuv(Ke), Kup(Ke),
421 Kvu(Ke), Kvv(Ke), Kvp(Ke),
422 Kpu(Ke), Kpv(Ke), Kpp(Ke);
433 std::vector<dof_id_type> dof_indices;
434 std::vector<dof_id_type> dof_indices_u;
435 std::vector<dof_id_type> dof_indices_v;
436 std::vector<dof_id_type> dof_indices_p;
437 std::vector<dof_id_type> dof_indices_alpha;
453 const Real theta = 1.;
463 for (
const auto & elem :
mesh.active_local_element_ptr_range())
473 dof_map.
dof_indices (elem, dof_indices_alpha, alpha_var);
475 const unsigned int n_dofs = dof_indices.size();
476 const unsigned int n_u_dofs = dof_indices_u.size();
477 const unsigned int n_v_dofs = dof_indices_v.size();
478 const unsigned int n_p_dofs = dof_indices_p.size();
484 fe_vel->reinit (elem);
485 fe_pres->reinit (elem);
493 Ke.
resize (n_dofs, n_dofs);
509 Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
510 Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
511 Kup.reposition (u_var*n_u_dofs, p_var*n_u_dofs, n_u_dofs, n_p_dofs);
513 Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
514 Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
515 Kvp.reposition (v_var*n_v_dofs, p_var*n_v_dofs, n_v_dofs, n_p_dofs);
517 Kpu.reposition (p_var*n_u_dofs, u_var*n_u_dofs, n_p_dofs, n_u_dofs);
518 Kpv.reposition (p_var*n_u_dofs, v_var*n_u_dofs, n_p_dofs, n_v_dofs);
519 Kpp.reposition (p_var*n_u_dofs, p_var*n_u_dofs, n_p_dofs, n_p_dofs);
522 Kp_alpha.reposition (p_var*n_u_dofs, p_var*n_u_dofs+n_p_dofs, n_p_dofs, 1);
523 Kalpha_p.reposition (p_var*n_u_dofs+n_p_dofs, p_var*n_u_dofs, 1, n_p_dofs);
526 Fu.reposition (u_var*n_u_dofs, n_u_dofs);
527 Fv.reposition (v_var*n_u_dofs, n_v_dofs);
528 Fp.reposition (p_var*n_u_dofs, n_p_dofs);
536 for (
unsigned int qp=0; qp<qrule.n_points(); qp++)
539 Number u = 0., u_old = 0.;
540 Number v = 0., v_old = 0.;
547 for (
unsigned int l=0; l<n_u_dofs; l++)
550 u_old += phi[l][qp]*navier_stokes_system.
old_solution (dof_indices_u[l]);
551 v_old += phi[l][qp]*navier_stokes_system.
old_solution (dof_indices_v[l]);
556 u += phi[l][qp]*navier_stokes_system.current_solution (dof_indices_u[l]);
557 v += phi[l][qp]*navier_stokes_system.current_solution (dof_indices_v[l]);
558 grad_u.
add_scaled (dphi[l][qp], navier_stokes_system.current_solution (dof_indices_u[l]));
559 grad_v.
add_scaled (dphi[l][qp], navier_stokes_system.current_solution (dof_indices_v[l]));
563 for (
unsigned int l=0; l<n_p_dofs; l++)
564 p_old += psi[l][qp]*navier_stokes_system.
old_solution (dof_indices_p[l]);
570 const Number u_x = grad_u(0);
571 const Number u_y = grad_u(1);
572 const Number v_x = grad_v(0);
573 const Number v_y = grad_v(1);
578 for (
unsigned int i=0; i<n_u_dofs; i++)
580 Fu(i) += JxW[qp]*(u_old*phi[i][qp] -
581 (1.-theta)*dt*(U_old*grad_u_old)*phi[i][qp] +
582 (1.-theta)*dt*p_old*dphi[i][qp](0) -
583 (1.-theta)*dt*nu*(grad_u_old*dphi[i][qp]) +
584 theta*dt*(U*grad_u)*phi[i][qp]);
587 Fv(i) += JxW[qp]*(v_old*phi[i][qp] -
588 (1.-theta)*dt*(U_old*grad_v_old)*phi[i][qp] +
589 (1.-theta)*dt*p_old*dphi[i][qp](1) -
590 (1.-theta)*dt*nu*(grad_v_old*dphi[i][qp]) +
591 theta*dt*(U*grad_v)*phi[i][qp]);
598 for (
unsigned int j=0; j<n_u_dofs; j++)
600 Kuu(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] +
601 theta*dt*nu*(dphi[i][qp]*dphi[j][qp]) +
602 theta*dt*(U*dphi[j][qp])*phi[i][qp] +
603 theta*dt*u_x*phi[i][qp]*phi[j][qp]);
605 Kuv(i,j) += JxW[qp]*theta*dt*u_y*phi[i][qp]*phi[j][qp];
607 Kvv(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] +
608 theta*dt*nu*(dphi[i][qp]*dphi[j][qp]) +
609 theta*dt*(U*dphi[j][qp])*phi[i][qp] +
610 theta*dt*v_y*phi[i][qp]*phi[j][qp]);
612 Kvu(i,j) += JxW[qp]*theta*dt*v_x*phi[i][qp]*phi[j][qp];
616 for (
unsigned int j=0; j<n_p_dofs; j++)
618 Kup(i,j) += JxW[qp]*(-theta*dt*psi[j][qp]*dphi[i][qp](0));
619 Kvp(i,j) += JxW[qp]*(-theta*dt*psi[j][qp]*dphi[i][qp](1));
627 for (
unsigned int i=0; i<n_p_dofs; i++)
629 Kp_alpha(i,0) += JxW[qp]*psi[i][qp];
630 Kalpha_p(0,i) += JxW[qp]*psi[i][qp];
631 for (
unsigned int j=0; j<n_u_dofs; j++)
633 Kpu(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](0);
634 Kpv(i,j) += JxW[qp]*psi[i][qp]*dphi[j][qp](1);
649 navier_stokes_system.rhs->add_vector (Fe, dof_indices);
654 navier_stokes_system.rhs->add(navier_stokes_system.rhs->size()-1, 10.);
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
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 &...)
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
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.