From: Derek Gaston <derek.gaston@in...>  20081031 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 offdiagonal. 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 