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