## [Libmesh-devel] Contraints for JFNK hanging nodes and Dirchlet BC's

 [Libmesh-devel] Contraints for JFNK hanging nodes and Dirchlet BC's From: Derek Gaston - 2008-10-31 18:57:03 ```So... I've been chatting with Roy and John, but I thought it was probably time to have an actual discussion about this. My immediate need is to apply Dirichlet BC constraints (strongly). I've been doing this by manually manipulating the element residual and Jacobian by finding if that element has a side on the boundary... but of course this won't work with tets (because they often have DOFs on a boundary and no face on that boundary)... and sometimes doesn't work for hexes and quads either. In order to get away from that I thought I would do a two step process: add all internal residual / jacobian contributions... global assemble the residual / jacobian... then loop over a set of sides or nodes and manually manipulate the entries. This isn't terrible... but not great. At the same time with JFNK we have a problem with the current way we deal with hanging node constraints in libMesh. For JFNK we need to fill the residual of a constrained DOF with something like u - ( A1 * Parent1Value + A2 * Parent2Value + etc.). Currently constrain_element_vector() doesn't do anything close to this (and if I'm reading it correctly) it actually zeroes out the constrained DOF's residual.... which means that nothing gets solved for, for that DOF. What we currently do is call enforce_constraints_exactly() immediately after the solve to fix up these values. Which is ok... but again not great. I'd rather just converge these values naturally along with the rest of the simulation. So... in light of these things I thought we might have a discussion about the constraint system. For the JFNK with hanging nodes problem I propose adding a new function: constrain_element_residual(). It will take everything constrain_element_vector() does now... plus a solution vector. Then it can write a residual like what I mentioned above into the correct spot in the residual. The Dirichlet BC problem gets more complicated. In my chatting with Roy it seems like it might be possible to add a new set of capabilities to the constraint system to handle Dirichlet BCs. I proposed an algorithm like: 1. Attach a compute_constraints() function to the system. 2. That compute_constraints() function would get called... it would loop over whatever it wanted to and add constraints to the dof_map.... along with the value to constrain the DOF to. 3. Assemble like usual... except that the constrain_element_matrix / vector stuff would automatically modify values to enforce the Dirichlet constraints. In a Newton update scheme you would supply the Dirichlet constraint value from compute_constraints() as u - value.... and that would get put into your RHS (residual) for that DOF. In the jacobian (if you have one) you would get a 1 on the diagonal and 0's in the off-diagonal. We could then extend this system to be a little easier for users. For instance, allowing them to pass a function of x,y,z that would compute the Dirichlet constraints or even specifying to pull the value out of a MeshData object or something. What do you guys think? Is there a better way? Derek ```

 [Libmesh-devel] Contraints for JFNK hanging nodes and Dirchlet BC's From: Derek Gaston - 2008-10-31 18:57:03 ```So... I've been chatting with Roy and John, but I thought it was probably time to have an actual discussion about this. My immediate need is to apply Dirichlet BC constraints (strongly). I've been doing this by manually manipulating the element residual and Jacobian by finding if that element has a side on the boundary... but of course this won't work with tets (because they often have DOFs on a boundary and no face on that boundary)... and sometimes doesn't work for hexes and quads either. In order to get away from that I thought I would do a two step process: add all internal residual / jacobian contributions... global assemble the residual / jacobian... then loop over a set of sides or nodes and manually manipulate the entries. This isn't terrible... but not great. At the same time with JFNK we have a problem with the current way we deal with hanging node constraints in libMesh. For JFNK we need to fill the residual of a constrained DOF with something like u - ( A1 * Parent1Value + A2 * Parent2Value + etc.). Currently constrain_element_vector() doesn't do anything close to this (and if I'm reading it correctly) it actually zeroes out the constrained DOF's residual.... which means that nothing gets solved for, for that DOF. What we currently do is call enforce_constraints_exactly() immediately after the solve to fix up these values. Which is ok... but again not great. I'd rather just converge these values naturally along with the rest of the simulation. So... in light of these things I thought we might have a discussion about the constraint system. For the JFNK with hanging nodes problem I propose adding a new function: constrain_element_residual(). It will take everything constrain_element_vector() does now... plus a solution vector. Then it can write a residual like what I mentioned above into the correct spot in the residual. The Dirichlet BC problem gets more complicated. In my chatting with Roy it seems like it might be possible to add a new set of capabilities to the constraint system to handle Dirichlet BCs. I proposed an algorithm like: 1. Attach a compute_constraints() function to the system. 2. That compute_constraints() function would get called... it would loop over whatever it wanted to and add constraints to the dof_map.... along with the value to constrain the DOF to. 3. Assemble like usual... except that the constrain_element_matrix / vector stuff would automatically modify values to enforce the Dirichlet constraints. In a Newton update scheme you would supply the Dirichlet constraint value from compute_constraints() as u - value.... and that would get put into your RHS (residual) for that DOF. In the jacobian (if you have one) you would get a 1 on the diagonal and 0's in the off-diagonal. We could then extend this system to be a little easier for users. For instance, allowing them to pass a function of x,y,z that would compute the Dirichlet constraints or even specifying to pull the value out of a MeshData object or something. What do you guys think? Is there a better way? Derek ```
 Re: [Libmesh-devel] Contraints for JFNK hanging nodes and Dirchlet BC's From: Benjamin Kirk - 2008-10-31 20:00:24 ```> At the same time with JFNK we have a problem with the current way we > deal with hanging node constraints in libMesh. For JFNK we need to > fill the residual of a constrained DOF with something like u - ( A1 * > Parent1Value + A2 * Parent2Value + etc.). Currently > constrain_element_vector() doesn't do anything close to this (and if > I'm reading it correctly) it actually zeroes out the constrained DOF's > residual.... which means that nothing gets solved for, for that DOF. In element residual computations, are you using DofMap::extract_local_vector()? http://libmesh.sourceforge.net/doxygen/classDofMap.php#cfbb9dbe71246f5d08607 7a3f4a5cced I use it to extract element-based Dofs, and they get properly constrained. Would the residual value then be the current value - the constraned value, as you mention? -Ben ```
 Re: [Libmesh-devel] Contraints for JFNK hanging nodes and Dirchlet BC's From: Benjamin Kirk - 2008-10-31 20:31:48 ``` >> At the same time with JFNK we have a problem with the current way we >> deal with hanging node constraints in libMesh. For JFNK we need to >> fill the residual of a constrained DOF with something like u - ( A1 * >> Parent1Value + A2 * Parent2Value + etc.). Currently >> constrain_element_vector() doesn't do anything close to this (and if >> I'm reading it correctly) it actually zeroes out the constrained DOF's >> residual.... which means that nothing gets solved for, for that DOF. > > In element residual computations, are you using > DofMap::extract_local_vector()? > http://libmesh.sourceforge.net/doxygen/classDofMap.php#cfbb9dbe71246f5d08607 > 7a3f4a5cced > > I use it to extract element-based Dofs, and they get properly constrained. > Would the residual value then be the current value - the constraned value, > as you mention? I guess what I am trying to say is than in ex19 rather than computing u & grad_u from soln passed in at the current nonlinear iterate, instead you would extract Ue from soln using the DofMap::extract_vector(), and that builds the constraint in. Then you compute u & grad_u from Ue. Then, for the constrained dofs, the residual value is simply (soln(dof_index[i]) - Ue(i))? That linear term falls through and gives you the desired action of the 1 in the Jacobian... -Ben ```
 Re: [Libmesh-devel] Contraints for JFNK hanging nodes and Dirchlet BC's From: Derek Gaston - 2008-10-31 20:38:08 ```On Oct 31, 2008, at 2:00 PM, Benjamin Kirk wrote: > I use it to extract element-based Dofs, and they get properly > constrained. > Would the residual value then be the current value - the constraned > value, > as you mention? I didn't know that method existed. It looks like it's pretty much achieving the same goal that enforce_constraints_exactly() does.... it at least means that you are calculating your next residual using the properly constrained value. Unfortunately, this doesn't help with calculating a constrained value. For instance, let's say that on solve #23 your solver decides it's converged... and libMesh writes out the solution. Do you have any idea what the value of your solution at the hanging nodes is? Unless you called enforce_constraints_exactly() before you wrote out the solution.... you don't. One idea would be for me to overwrite my own residual before I add it into the global residual.... I guess I could call is_constrained() on every node on the element and if it is then I could overwrite what residual is there with u - constrained_value before I add it in..... this is essentially the work I was thinking constrain_element_residual() would do. Derek ```