From: Roy Stogner <roystgnr@ic...>  20060119 15:12:30

On Thu, 19 Jan 2006, Shengli Xu wrote: > 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? Yes. > 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 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_indices); This function is to handle hanging degrees of freedom on adaptively refined meshes. With first order elements, for example, the degree of freedom at a "hanging node" can't take just any value  it has to be the average of the values of the DoFs at the two neighboring nodes. This shouldn't be a concern for different boundary condition specifications; hanging nodes are on sides between differently refined elements and so by definition aren't on the boundary. > I checked that the global stiffness matrix and right load vector are wright. I'm assuming thie means "right, before boundary conditions are added". If you get the right matrix and force vector but the wrong solution, then you've got a solver problem. > 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. This code will break horribly on nonLagrange nonfirst order elements, but should be fine for the elements you're using, assuming your geometric elements are also first order and once you find whatever the bug is. This looks basically like what we do in examples 6, 7, and 8. Have you printed out postboundarycondition matrices to examine? It should be easy to make sure the right number of DoFs are getting attached penalty terms, and then not too much harder to make sure the right DoFs are getting them.  Roy Stogner 