156 libmesh_assert_equal_to (system_name,
"Stokes");
182 std::unique_ptr<FEBase> fe_vel (FEBase::build(
dim, fe_vel_type));
186 std::unique_ptr<FEBase> fe_pres (FEBase::build(
dim, fe_pres_type));
193 fe_vel->attach_quadrature_rule (&qrule);
194 fe_pres->attach_quadrature_rule (&qrule);
200 const std::vector<Real> & JxW = fe_vel->get_JxW();
204 const std::vector<std::vector<RealGradient>> & dphi = fe_vel->get_dphi();
208 const std::vector<std::vector<Real>> & psi = fe_pres->get_phi();
224 Kuu(Ke), Kuv(Ke), Kup(Ke),
225 Kvu(Ke), Kvv(Ke), Kvp(Ke),
226 Kpu(Ke), Kpv(Ke), Kpp(Ke);
236 std::vector<dof_id_type> dof_indices;
237 std::vector<dof_id_type> dof_indices_u;
238 std::vector<dof_id_type> dof_indices_v;
239 std::vector<dof_id_type> dof_indices_p;
249 for (
const auto & elem :
mesh.active_local_element_ptr_range())
260 const unsigned int n_dofs = dof_indices.size();
261 const unsigned int n_u_dofs = dof_indices_u.size();
262 const unsigned int n_v_dofs = dof_indices_v.size();
263 const unsigned int n_p_dofs = dof_indices_p.size();
269 fe_vel->reinit (elem);
270 fe_pres->reinit (elem);
278 Ke.
resize (n_dofs, n_dofs);
294 Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs);
295 Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs);
296 Kup.reposition (u_var*n_u_dofs, p_var*n_u_dofs, n_u_dofs, n_p_dofs);
298 Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs);
299 Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs);
300 Kvp.reposition (v_var*n_v_dofs, p_var*n_v_dofs, n_v_dofs, n_p_dofs);
302 Kpu.reposition (p_var*n_u_dofs, u_var*n_u_dofs, n_p_dofs, n_u_dofs);
303 Kpv.reposition (p_var*n_u_dofs, v_var*n_u_dofs, n_p_dofs, n_v_dofs);
304 Kpp.reposition (p_var*n_u_dofs, p_var*n_u_dofs, n_p_dofs, n_p_dofs);
306 Fu.reposition (u_var*n_u_dofs, n_u_dofs);
307 Fv.reposition (v_var*n_u_dofs, n_v_dofs);
308 Fp.reposition (p_var*n_u_dofs, n_p_dofs);
311 for (
unsigned int qp=0; qp<qrule.n_points(); qp++)
315 for (
unsigned int i=0; i<n_u_dofs; i++)
316 for (
unsigned int j=0; j<n_u_dofs; j++)
317 Kuu(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
320 for (
unsigned int i=0; i<n_u_dofs; i++)
321 for (
unsigned int j=0; j<n_p_dofs; j++)
322 Kup(i,j) += -JxW[qp]*psi[j][qp]*dphi[i][qp](0);
327 for (
unsigned int i=0; i<n_v_dofs; i++)
328 for (
unsigned int j=0; j<n_v_dofs; j++)
329 Kvv(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
332 for (
unsigned int i=0; i<n_v_dofs; i++)
333 for (
unsigned int j=0; j<n_p_dofs; j++)
334 Kvp(i,j) += -JxW[qp]*psi[j][qp]*dphi[i][qp](1);
339 for (
unsigned int i=0; i<n_p_dofs; i++)
340 for (
unsigned int j=0; j<n_u_dofs; j++)
341 Kpu(i,j) += -JxW[qp]*psi[i][qp]*dphi[j][qp](0);
344 for (
unsigned int i=0; i<n_p_dofs; i++)
345 for (
unsigned int j=0; j<n_v_dofs; j++)
346 Kpv(i,j) += -JxW[qp]*psi[i][qp]*dphi[j][qp](1);
365 for (
auto s : elem->side_index_range())
366 if (elem->neighbor_ptr(s) ==
nullptr)
368 const Elem & side = side_builder(*elem, s);
380 const Real penalty = 1.e10;
385 const Real u_value = (yf > .99) ? 1. : 0.;
388 const Real v_value = 0.;
393 for (
auto n : elem->node_index_range())
394 if (elem->node_id(n) == side.
node_id(ns))
401 Fu(n) += penalty*u_value;
402 Fv(n) += penalty*v_value;
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
void constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix and vector.
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and linear solvers ...
void resize(const unsigned int n)
Resize the vector.
virtual void add_vector(const T *v, const std::vector< numeric_index_type > &dof_indices)
Computes , where v is a pointer and each dof_indices[i] specifies where to add value v[i]...
This is the base class from which all geometric element types are derived.
NumericVector< Number > * rhs
The system matrix.
Order default_quadrature_order() const
const T_sys & get_system(std::string_view name) const
This is the MeshBase class.
unsigned int variable_number(std::string_view var) const
Defines a dense subvector for use in finite element computations.
This class handles the numbering of degrees of freedom on a mesh.
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.
Defines a dense submatrix for use in Finite Element-type computations.
Helper for building element sides that minimizes the construction of new elements.
const FEType & variable_type(const unsigned int i) 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
IntRange< unsigned short > node_index_range() const
const DofMap & get_dof_map() const
dof_id_type node_id(const unsigned int i) const
const Point & point(const unsigned int i) const
const SparseMatrix< Number > & get_system_matrix() const