Index: ex13.C =================================================================== --- ex13.C (revision 3410) +++ ex13.C (working copy) @@ -120,6 +120,10 @@ // providing an LBB-stable pressure-velocity pair. system.add_variable ("p", FIRST); + // Add a scalar Lagrange multiplier to constrain the + // pressure to have zero mean. + system.add_variable ("alpha", FIRST, SCALAR); + // Give the system a pointer to the matrix assembly // function. system.attach_assemble_function (assemble_stokes); @@ -295,6 +299,7 @@ const unsigned int u_var = navier_stokes_system.variable_number ("u"); const unsigned int v_var = navier_stokes_system.variable_number ("v"); const unsigned int p_var = navier_stokes_system.variable_number ("p"); + const unsigned int alpha_var = navier_stokes_system.variable_number ("alpha"); // Get the Finite Element type for "u". Note this will be // the same as the type for "v". @@ -356,6 +361,7 @@ Kuu(Ke), Kuv(Ke), Kup(Ke), Kvu(Ke), Kvv(Ke), Kvp(Ke), Kpu(Ke), Kpv(Ke), Kpp(Ke); + DenseSubMatrix Kalpha_p(Ke), Kp_alpha(Ke); DenseSubVector Fu(Fe), @@ -369,6 +375,7 @@ std::vector dof_indices_u; std::vector dof_indices_v; std::vector dof_indices_p; + std::vector dof_indices_alpha; // Find out what the timestep size parameter is from the system, and // the value of theta for the theta method. We use implicit Euler (theta=1) @@ -406,6 +413,7 @@ dof_map.dof_indices (elem, dof_indices_u, u_var); dof_map.dof_indices (elem, dof_indices_v, v_var); dof_map.dof_indices (elem, dof_indices_p, p_var); + dof_map.dof_indices (elem, dof_indices_alpha, alpha_var); const unsigned int n_dofs = dof_indices.size(); const unsigned int n_u_dofs = dof_indices_u.size(); @@ -453,6 +461,11 @@ Kpv.reposition (p_var*n_u_dofs, v_var*n_u_dofs, n_p_dofs, n_v_dofs); Kpp.reposition (p_var*n_u_dofs, p_var*n_u_dofs, n_p_dofs, n_p_dofs); + // Also, add a row and a column to constrain the pressure + Kp_alpha.reposition (p_var*n_u_dofs, p_var*n_u_dofs+n_p_dofs, n_p_dofs, 1); + Kalpha_p.reposition (p_var*n_u_dofs+n_p_dofs, p_var*n_u_dofs, 1, n_p_dofs); + + Fu.reposition (u_var*n_u_dofs, n_u_dofs); Fv.reposition (v_var*n_u_dofs, n_v_dofs); Fp.reposition (p_var*n_u_dofs, n_p_dofs); @@ -558,6 +571,8 @@ // negative one. Here we do not. for (unsigned int i=0; ineighbor(side) == NULL) - // Pin the pressure to zero at global node number "pressure_node". - // This effectively removes the non-trivial null space of constant - // pressure solutions. - const bool pin_pressure = true; - if (pin_pressure) - { - const unsigned int pressure_node = 0; - const Real p_value = 0.0; - for (unsigned int c=0; cn_nodes(); c++) - if (elem->node(c) == pressure_node) - { - Kpp(c,c) += penalty; - Fp(c) += penalty*p_value; - } - } +// // Pin the pressure to zero at global node number "pressure_node". +// // This effectively removes the non-trivial null space of constant +// // pressure solutions. +// const bool pin_pressure = true; +// if (pin_pressure) +// { +// const unsigned int pressure_node = 0; +// const Real p_value = 0.0; +// for (unsigned int c=0; cn_nodes(); c++) +// if (elem->node(c) == pressure_node) +// { +// Kpp(c,c) += penalty; +// Fp(c) += penalty*p_value; +// } +// } } // end boundary condition section - + // If this assembly program were to be used on an adaptive mesh, // we would have to apply any hanging node constraint equations dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices); @@ -654,7 +669,10 @@ navier_stokes_system.matrix->add_matrix (Ke, dof_indices); navier_stokes_system.rhs->add_vector (Fe, dof_indices); } // end of element loop - + + // We can set the mean of the pressure by setting Falpha + navier_stokes_system.rhs->add(navier_stokes_system.rhs->size()-1,10.); + // That's it. return; }