In example 11, the boundary condition is addressed in the loop of assembling the stiffness and right hand force. Can the boundary condition be addressed after  the global  stiffness  matrix and the global right load vector have assembled?
I think they are equal. But what does the following code mean?

We have now built the element matrix and RHS vector in terms of the element degrees of freedom. However, it is possible that some of the element DOFs are constrained to enforce solution continuity, i.e. they are not really "free". We need to constrain those DOFs in terms of non-constrained DOFs to ensure a continuous solution. The \p DofMap::constrain_element_matrix_and_vector() method does just that.
              dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);

In my program, the global stiffness matrix and the right load vector are assembled first. Then the boundar conditions
of u,v,p are addresed through the loop of node. But I can not get the correct result. My problem is to solve 2D stokes equation using
the least square finite element. Some information about my program is :
// add variable "u","v","p","oz"
system.add_variable("u",FIRST);
system.add_variable("v",FIRST);
system.add_variable("p",FIRST);
system.add_variable("oz",FIRST); // oz is the vorticity

I checked that the global stiffness matrix and right load vector are wright.

After assembled the stiffness matrix and right load vector, the boundary condition is addressed as: (the domain is 1X1,
the boundary condition is u=0,v=0 at the four sides of domain, p=0 at node (0.5, 0).
the code is:
/////////////////////////////////
for(unsigned int n_cnt=0;n_cnt<n_nodes;n_cnt++){
const Node &curr_node = mesh.node(n_cnt);

// set boundary on y=0
const Real y_bound = 0.;
if(fabs(curr_node(1)-y_bound) < TOLERANCE){
unsigned int dn = curr_node.dof_number(0,0,0);
const Real penalty = 1.e10;
const Real u_value = 0.0;
unsigned int dny = curr_node.dof_number(0,1,0);
const Real v_value = 0.0;
if(fabs(curr_node(0)-0.5) < TOLERANCE){
unsigned int dnp = curr_node.dof_number(0,2,0);
const Real p_value = 0.0;
rhs.add(dnp,p_value*penalty);
matrix.add(dnp,dnp,penalty);
}

rhs.add(dn,u_value*penalty);
matrix.add(dn,dn,penalty);

rhs.add(dny,v_value*penalty);
matrix.add(dny,dny,penalty);
}

// set boundary on y=1.
const Real y_top = 1.;
if(fabs(curr_node(1)-y_top) < TOLERANCE){
unsigned int dn = curr_node.dof_number(0,0,0);
const Real u_value = 0.0;
unsigned int dny = curr_node.dof_number(0,1,0);
const Real penalty = 1.e10;
const Real v_value = 0.0;

rhs.add(dn,u_value*penalty);
matrix.add(dn,dn,penalty);

rhs.add(dny,v_value*penalty);
matrix.add(dny,dny,penalty);
}
// set boundary on x=0.
const Real x_left = 0.;
if(fabs(curr_node(0)-x_left) < TOLERANCE && (fabs(curr_node(1)-y_bound)>TOLERANCE) && (fabs(curr_node(1)-y_top)>TOLERANCE) ){
unsigned int dn = curr_node.dof_number(0,0,0);
const Real penalty = 1.e10;
const Real u_value = 0.0;
unsigned int dny =curr_node.dof_number(0,1,0);
const Real v_value = 0.0;

rhs.add(dn,u_value*penalty);
matrix.add(dn,dn,penalty);

rhs.add(dny,v_value*penalty);
matrix.add(dny,dny,penalty);
}
// set boundary on x=1.
const Real x_right = 1.;
if(fabs(curr_node(0)-x_right) < TOLERANCE &&(fabs(curr_node(1)-y_bound)>TOLERANCE) && (fabs(curr_node(1)-y_top)>TOLERANCE)){
unsigned int dn = curr_node.dof_number(0,0,0);
const Real penalty = 1.e10;
const Real u_value = 0.0;
unsigned int dny = curr_node.dof_number(0,1,0);
const Real v_value = 0.0;

rhs.add(dn,u_value*penalty);
matrix.add(dn,dn,penalty);

rhs.add(dny,v_value*penalty);
matrix.add(dny,dny,penalty);
}
////////////////////////////
I want to know is it right to address boundary condition.

ShengliXu