## Re: [Libmesh-users] Questions about example 11

 Re: [Libmesh-users] Questions about example 11 From: Roy Stogner - 2006-01-19 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 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); 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 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 non-Lagrange non-first 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 post-boundary-condition 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 ```

 [Libmesh-users] Questions about example 11 From: Shengli Xu - 2006-01-19 13:46:56 Attachments: Message as HTML ```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 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_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_cntTOLERANCE) && (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 ```
 Re: [Libmesh-users] Questions about example 11 From: Roy Stogner - 2006-01-19 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 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); 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 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 non-Lagrange non-first 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 post-boundary-condition 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 ```