From: Shengli Xu <shengli.xu.xu@gm...>  20060119 13:46:56

In example 11, the boundary condition is addressed in the loop of assemblin= g 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 elemen= t 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 nonconstrained 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_indi= ces); 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=3D0,v=3D0 at the four sides of domain, p=3D0 at node (0.5, 0). the code is: ///////////////////////////////// for(unsigned int n_cnt=3D0;n_cnt<n_nodes;n_cnt++){ const Node &curr_node =3D mesh.node(n_cnt); // set boundary on y=3D0 const Real y_bound =3D 0.; if(fabs(curr_node(1)y_bound) < TOLERANCE){ unsigned int dn =3D curr_node.dof_number(0,0,0); const Real penalty =3D 1.e10; const Real u_value =3D 0.0; unsigned int dny =3D curr_node.dof_number(0,1,0); const Real v_value =3D 0.0; if(fabs(curr_node(0)0.5) < TOLERANCE){ unsigned int dnp =3D curr_node.dof_number(0,2,0); const Real p_value =3D 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=3D1. const Real y_top =3D 1.; if(fabs(curr_node(1)y_top) < TOLERANCE){ unsigned int dn =3D curr_node.dof_number(0,0,0); const Real u_value =3D 0.0; unsigned int dny =3D curr_node.dof_number(0,1,0); const Real penalty =3D 1.e10; const Real v_value =3D 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=3D0. const Real x_left =3D 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 =3D curr_node.dof_number(0,0,0); const Real penalty =3D 1.e10; const Real u_value =3D 0.0; unsigned int dny =3Dcurr_node.dof_number(0,1,0); const Real v_value =3D 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=3D1. const Real x_right =3D 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 =3D curr_node.dof_number(0,0,0); const Real penalty =3D 1.e10; const Real u_value =3D 0.0; unsigned int dny =3D curr_node.dof_number(0,1,0); const Real v_value =3D 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 